Topics in PhD Level Econometrics with Streamlit
import streamlit as st
from st_pages import Page, add_page_title, show_pages_from_config

import src.scripts.plot_themes
import src.scripts.utils as utl

st.set_page_config(
    page_title="PhD Econometrics",
    page_icon="📈",
    layout="wide",
)

utl.local_css("src/styles/styles_home.css")
utl.external_css(
    "https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.15.3/css/all.min.css"
)

show_pages_from_config()

s1, c1, c2 = utl.wide_col()

# my LinkedIn, GitHub, and email
linkedin_url = "https://www.linkedin.com/in/justinas-grigaitis/"
github_url = "https://github.com/justgri"
email_url = "mailto:justinas.grigaitis@econ.uzh.ch"

# Intro
with c1:
    # Title iterations
    # Econometrics for PhD Students
    # PhD Econometrics for Everyone
    # PhD Econometrics for All
    # PhD for All: Econometrics

    st.title("PhD for All: Econometrics")

    # Subheader iterations
    # Intuitive visuals. Rigorous theory. Reproducible code.
    # Intuitive visuals. Rigorous theory. Simple code.
    # Intuitive visuals. Rigorous theory. Easy code.
    # Interactive visuals. Rigorous theory. Simple code.

    st.markdown(
        '<span style="font-size: 28px; display: block; margin-bottom: 5px;">*Interactive visuals. Rigorous theory. Simple code.*</span>',
        unsafe_allow_html=True,
    )

    st.markdown(
        "<hr style='margin-top: 0; margin-bottom: 5px;'>",
        unsafe_allow_html=True,
    )

    st.markdown(
        r"""Learning and helping others learn along the way.<br>
            Explaining PhD concepts intuitively.<br>
            Procrastinating productively.""",
        unsafe_allow_html=True,
    )
    st.markdown(
        r"""Targeted at **grad students**, but useful for **professionals** and **undergrads** alike.""",
        unsafe_allow_html=True,
    )
    st.markdown(
        f"""
        Please send me feedback:<br>
    <a href="{linkedin_url}" target="_blank">
        <i class="fab fa-linkedin fa-lg"></i>
    </a>
    <a href="{email_url}" target="_blank">
        <i class="fas fa-envelope fa-lg"></i>
    </a>
    <a href="{github_url}" target="_blank">
        <i class="fab fa-github fa-lg"></i>
    </a>
    """,
        unsafe_allow_html=True,
    )

    st.markdown(
        r"""<u>**Disclaimer:**</u> <br>
        This website does not represent the official curriculum taught at my university. <br>
        My goal is to cover fewer topics in greater depth rather than scratch the surface of many. <br>
        All mistakes are my own.
        """,
        # Main difference is matrix algebra and proving everything along the way, which might not always be included here.
        # Hopefully it will give insights to both PhD students, undergrads, and others.
        unsafe_allow_html=True,
    )


_, note_col, _ = st.columns((0.02, 1, 1.5))
with note_col:
    with st.expander("For mobile users:"):
        st.write("Sidebar leads to other pages within this app.")
        st.image("src/images/intro_mobile_sidebar.png", width=200)

s1, c2, s2 = utl.wide_col()

# Preliminary ToC
with c2:
    st.markdown(
        "<h3 style='text-align: center'>Table of Contents</h3>",
        unsafe_allow_html=True,
    )

    with st.expander("Click to expand"):
        # Page links - potentially hrefs with st.experimental_set_query_params()
        path_tools = "https://phd-econometrics.streamlit.app/Must-know"
        path_ols = "https://phd-econometrics.streamlit.app/Linear%20Regression"
        path_fit = "https://phd-econometrics.streamlit.app/Fit%20Measures"

        st.markdown(
            r"""
            Hyperlinks lead to the corresponding pages on this website.""",
            unsafe_allow_html=True,
        )

        st.markdown(
            f"""[**Top 10 theory tools**]({path_tools})"""
            + r""" **that everyone should know according to Jeffrey Wooldridge**
        """,
            unsafe_allow_html=True,
        )

        st.markdown(
            rf"""
        <div class="numbered-header">
            <b>Chapter 1: Finite-Sample Properties of OLS</b><br>
        </div>
        
        <div class="numbered">
            1. <a href="{path_ols}" target="_blank">OLS estimation</a> (Greene Ch 2 - 3.2)<br>
            2. <a href="{path_fit}" target="_blank">Fit measures</a> (Greene Ch 3.5 and Ch 5.8) <br>
            3. Hypothesis testing (Greene Ch 5) <br>
            4. Data problems (OVB, measurement error, missing data - Greene Ch. 4.9)<br>
            5. Functional form (Greene Ch 6.5)<br>
            6. PCA (Greene Ch 4.9)
        </div>

        <br>

        <div class="numbered-header">
            <b>Chapter 2: Large-Sample Theory</b><br>
        </div>

        <div class="numbered">
            7. Limit theorems, "delta-method" <br>
            8. Law of large numbers <br>
            9. Large sample OLS properties <br>
        </div>

        <br>

        <div class="numbered-header">
            <b>Chapter 3: Generalized Method of Moments</b><br>
        </div>

        <div class="numbered">
            10. Endogeneity <br>
            11. Instrumental Variables <br>
            12. IV-related tests <br>
            13. GMM properties
        </div>

        <br>
        
        <div class="numbered-header">
            <b>Chapter 4: Multiple-Equation GMM</b>
        </div>

        <div class="numbered">
            14. TBD <br>
            15. TBD <br>
        </div>

        <br>

        
        
        Next semester - Panel Data, Time Series, Cointegration, and MLE. <br>
        Bonus data science - causal ML, deep learning, classification algorithms, etc.

        Section headers follow Hayahshi *Econometrics* (1st ed.).<br>
        Chapters from Greene *Eceonometric Analysis* (8th ed.) given in the parentheses.<br>
        Subsections are likely to change depending on which topics I find most interesting or challenging.<br>

        """,
            unsafe_allow_html=True,
        )


with c2:
    st.markdown(
        "<h3 style='text-align: center'>What is Econometrics?</h3>",
        unsafe_allow_html=True,
    )

    econometrica_public = "https://www.sv.uio.no/econ/om/tall-og-fakta/nobelprisvinnere/ragnar-frisch/published-scientific-work/rf-published-scientific-works/rf1933c.pdf"
    econometrica_jstor = "https://www.jstor.org/stable/i332704"

    with st.expander("Click to expand"):
        st.markdown(
            r"""
            "\[E\]conometrics is by no means the same as economic statistics.
            Nor is it identical with what we call general economic theory,
            although a considerable portion of this theory has a definitely quantitative character.
            Nor should econometrics be taken as synonomous with the application of mathematics to economics.
            Experience has shown that each of these three view-points, that of statistics, economic theory, and mathematics, is a necessary, but not by itself a sufficient,
            condition for a real understanding of the quan- titative relations in modern economic life.
            **It is the unification of all three that is powerful.
            And it is this *unification* that constitutes econometrics.**" (emphasis added)
            """,
            unsafe_allow_html=True,
        )

        st.markdown(
            f"""Frisch, Ragnar. “Editor’s Note.” Econometrica, vol. 1, no. 1, 1933, pp. 1–4. JSTOR, http://www.jstor.org/stable/1912224. Accessed 1 Oct. 2023. (free access [link]({econometrica_public}))""",
            unsafe_allow_html=True,
        )


# Textbooks
with c2:
    st.markdown(
        "<h3 style='text-align: center'>Reference Textbooks</h3>",
        unsafe_allow_html=True,
    )

    c2_1, s2_1, c2_2 = st.columns((1, 0.05, 1))

    with c2_1:
        # st.image("src/images/intro_stock_watson.jpg", width=350)
        st.image("src/images/intro_hayashi.jpg", width=300)

    with c2_2:
        st.image("src/images/intro_greene.jpg", width=300)


# Other references
with c2:
    st.markdown(
        "<h3 style='text-align: center'>Other References</h3>",
        unsafe_allow_html=True,
    )

    st.link_button(
        "NYU Lecture Notes by William H. Greene",
        "https://pages.stern.nyu.edu/~wgreene/Econometrics/Notes.htm",
        type="secondary",
    )

    st.link_button(
        "MIT Lecture Notes by Victor Chernozhukov",
        "https://ocw.mit.edu/courses/14-382-econometrics-spring-2017/pages/lecture-notes/",
        type="secondary",
    )

    st.link_button(
        "An Introduction to Statistical Learning with R/Python (free textbooks)",
        "https://www.statlearning.com/",
        type="secondary",
    )

    st.link_button(
        "TA notes by Matteo Courthoud",
        "https://matteocourthoud.github.io/course/metrics/",
        type="secondary",
    )

    st.link_button(
        "Graduate Econometrics YouTube playlist by Ben Lambert",
        "https://www.youtube.com/watch?v=GMVh02WGhoc&list=PLwJRxp3blEvaxmHgI2iOzNP6KGLSyd4dz",
        type="secondary",
    )

    st.link_button(
        "Intro to Econometrics with R by Professors at University of Duisburg-Essen",
        "https://www.econometrics-with-r.org/",
        type="secondary",
    )

    st.link_button(
        "Intro to Econometrics with R by Professors at SciencesPo - more focus on practical examples and intuition",
        "https://scpoecon.github.io/ScPoEconometrics/index.html",
        type="secondary",
    )

intro.py
justgri/phd_econometrics
import numpy as np
import src.scripts.plot_themes as thm
import src.scripts.utils as utl
import streamlit as st
from matplotlib import pyplot as plt
from st_pages import add_page_title

### PAGE CONFIGS ###

st.set_page_config(
    page_title="PhD Econometrics - Top Tools",
    page_icon="📈",
    layout="wide",
)

utl.local_css("src/styles/styles_pages.css")

random_seed = 0
s1, c1, s2 = utl.wide_col()

### PAGE START ###
# Dashboard header
with c1:
    st.title("Top 10 Theory Tools for Econometrics")
    st.divider()
    st.markdown(r"""<h3>According to Jeffrey Wooldridge</h3>""", unsafe_allow_html=True)


s1, c2, s2 = utl.wide_col()

#### EXAMPLE CHART FOR LINEARITY OF EXPECTATION AND VARIANCE
if "random_seed" not in st.session_state:
    st.session_state.random_seed = 0


def render_lin_exp():
    header_col = st.container()

    def generate_data(a, b, expected_value):
        np.random.seed(st.session_state.random_seed)

        # Generate the first deviation (always positive)
        deviation1 = np.random.uniform(0.3, 1)
        deviation2 = np.random.uniform(deviation1 + 0.5, deviation1 + 1.5)
        deviations = np.array([deviation1, deviation2])

        # Calculate mirror deviations
        mirror_deviations = -deviations

        # Combine and adjust for the expected value
        X = np.concatenate([deviations, mirror_deviations]) + expected_value
        Y = a + b * X

        return X, Y

    def format_two_decimal(x):
        rounded_x = round(x, 2)
        return f"{rounded_x:.2f}" if abs(rounded_x) > 1e-6 else "0.00"

    def plot_expectation(X, Y, a, b):
        fig, ax = plt.subplots(figsize=(4, 4))
        plt.subplots_adjust(left=0)  # remove margin
        plt.xlabel("X", fontweight="bold")
        plt.ylabel("Y", fontweight="bold")
        y_lb, y_ub = -8, 8
        x_lb, x_ub = -4, 4

        ax.yaxis.set_label_coords(-0.08, 0.5)

        ax.scatter(X, Y, color=thm.set1_blue)

        # Regression Line (Deterministic Line based on the given a and b)
        x_values = [min(X), max(X)]  # first and last X
        y_values = a + b * np.array(x_values)
        ax.plot(x_values, y_values, color=thm.set1_orange, label=r"$Y = a + bX$")

        # X BAR line
        E_X = np.mean(X)
        ax.axvline(E_X, color=thm.set1_blue, linestyle="-", linewidth=2)

        ### X BAR text
        x_offset = -0.18
        ax.text(
            E_X + x_offset,
            7.0,
            r"$\bar{X}=$" + format_two_decimal(E_X),
            color=thm.set1_blue,
            verticalalignment="center",
            horizontalalignment="right",
        )

        # Y BAR line
        E_Y = np.mean(Y)
        ax.axhline(E_Y, color=thm.set1_green, linestyle="-", linewidth=2)

        # Y BAR text
        y_offset = 0.5
        ax.text(
            -3.9,
            E_Y + y_offset,
            r"$\bar{Y}=$" + format_two_decimal(E_Y),
            color=thm.set1_green,
            verticalalignment="center",
        )

        # Dashed lines from x_bar to each data point
        has_lab_x_res = False
        var_X = np.var(X, ddof=0)

        for xi, yi in zip(X, Y):
            if not has_lab_x_res:
                lab_x_res = r"$Var(X)=$" + format_two_decimal(var_X)
                has_lab_x_res = True
            else:
                lab_x_res = None  # No label for other lines
            ax.plot(
                [E_X, xi],
                [yi, yi],
                color=thm.set1_blue,
                linestyle="--",
                linewidth=1,
                label=lab_x_res,
            )

        # Dashed lines from y_bar to each data point
        has_lab_y_res = False
        var_Y = np.var(Y, ddof=0)

        for xi, yi in zip(X, Y):
            if not has_lab_y_res:
                lab_y_res = r"$Var(Y)=$" + format_two_decimal(var_Y)
                has_lab_y_res = True
            else:
                lab_y_res = None  # No label for other lines

            ax.plot(
                [xi, xi],
                [E_Y, yi],
                color=thm.set1_green,
                linestyle="--",
                linewidth=1,
                label=lab_y_res,
            )

        # Update the y-axis limits after everything has been plotted
        ax.set_ylim(y_lb, y_ub)
        ax.set_xlim([x_lb, x_ub])
        ax.legend(loc="lower left", frameon=True, fontsize="small")

        return fig

    slider_col, chart_col = st.columns((0.75, 1))

    with slider_col:
        st.write("**Choose the following parameters:**")

        E_X_slider = st.slider(
            r"Expectation $E[X]$",
            min_value=-1.0,
            max_value=1.0,
            value=0.5,
            step=0.1,
        )

        a_slider = st.slider(
            r"Intercept $a$",
            min_value=-2.0,
            max_value=2.0,
            value=0.0,
            step=0.1,
        )
        b_slider = st.slider(
            r"Slope $b$", min_value=-2.0, max_value=2.0, value=1.0, step=0.1
        )

        X, Y = generate_data(a_slider, b_slider, E_X_slider)

        # Button to resample X
        if st.button("Resample X", type="primary"):
            st.session_state.random_seed += 1  # Increment seed by 1
            X, Y = generate_data(a_slider, b_slider, E_X_slider)

    with chart_col:
        st.pyplot(plot_expectation(X, Y, a_slider, b_slider))

    with header_col:
        if b_slider < 0:
            b_str = f"{b_slider:.2f}"
        else:
            b_str = f"+ {b_slider:.2f}"

        var_X_calc = np.var(X, ddof=0)  # ddof=0 ensures we get the population variance
        var_Y_calc = b_slider**2 * var_X_calc

        #         st.markdown(r"""$E[Y] = E[a+bX] = a + bE[X]$""", unsafe_allow_html=True)
        # st.markdown(r"""$Var(Y) = Var(a+bX) = b^2Var(X)$""", unsafe_allow_html=True)

        exp_eqn = (
            rf"$E[Y] = E[a+bX] = a + bE[X] = "
            + r"a + b \frac{\sum_{i=1}^N X_i}{N} ="
            + rf"{a_slider:.2f} {b_str} \times {E_X_slider:.2f} = {a_slider + b_slider * E_X_slider:.2f}$"
        )
        var_eqn = (
            rf"""$Var(Y) = Var(a+bX) = b^2Var(X) = ({b_slider:.2f})^2 """
            + r"\frac{\sum_{i=1}^N(X_i - E[X])^2}{N}"
            + rf"= {b_slider**2:.2f} \times {var_X_calc:.2f} = {var_Y_calc:.2f}$"
            ""
        )

        st.markdown(exp_eqn, unsafe_allow_html=True)
        st.markdown(var_eqn, unsafe_allow_html=True)


# Wooldridge Top 10
with c2:
    jmw_source = "https://docs.google.com/document/d/1GQ3qlD0_cNQdkh_iT-dBDXwNVIXk27WCZtmQ3Sgxz1o/mobilebasic#id.f68ahvvbxdaj"
    jmw_twitter = "https://twitter.com/jmwooldridge"
    tr_source = "https://bookdown.org/ts_robinson1994/10EconometricTheorems/"
    st.write(
        f"Source: Jeffrey Wooldridge [@jmwooldridge]({jmw_twitter}). For his tweets organized by topic check this - [link]({jmw_source})."
    )
    st.write(
        f"Thomas Robinson from LSE elaborate on each of them in this nice e-book - [link]({tr_source})."
    )

    # st.markdown(
    #     "<hr style='margin-top: 0; margin-bottom: 5px;'>",
    #     unsafe_allow_html=True,
    # )

    st.markdown(
        r"""
    <h5>1. Law of Iterated Expectations, Law of Total Variance </h5>
    
    $E[x] = E[E[x|y]]$ <br>
    $var[Y] = E[var[Y|X]] + var(E[Y|X])$

    <h5>2. Linearity of Expectations, Variance of a Sum</h5>

    $E[aX+bY+c] = aE[X] + bE[Y] + c$ <br>
    $var[aX+bY+c] = a^2 var(X) + b^2 var(Y) + 2ab Cov(X,Y)$
    """,
        unsafe_allow_html=True,
    )

st.markdown(
    r"""<h4 style="text-align:center;">Visualize a Simple Example</h4>""",
    unsafe_allow_html=True,
)

st.markdown(
    r"""
    Consider a random variable $X$ and a random variable $Y = a + bX$, i.e., a linear function of $X$.<br>
    Move the sliders and observe how $Var(Y)$ only depends on $b$, but not on $E[X]$ or $a$,
    whereas $E[Y]$ changes linearly based on the combination of all three parameters.
    Also note that $E[Y]$ does not depend on $b$ if $E[X] = 0$ and vice-versa.""",
    unsafe_allow_html=True,
)

render_lin_exp()

s1, c3, s2 = utl.wide_col()

with c3:
    st.markdown(
        r"""
    <h5>3. Jensen's Inequality, Chebyshev's Inequality</h5>

    Jensen:
    $\varphi(\operatorname{E}[X]) \leq \operatorname{E} \left[\varphi(X)\right]$, where $X$ is a RV and $\varphi$ is convex <br>
    Chebyshev:
    $\Pr(|X-\mu|\geq k\sigma) \leq \frac{1}{k^2}$, where $X$ is a RV, $σ^2$ is finite, and $k>0$

    <h5>4. Linear Projection and its Properties</h5>

    $\mathbf{P} \equiv \mathbf{X}(\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'$ and $\mathbf{M} \equiv \mathbf{I}_n - \mathbf{P}$ <br> 
    $\mathbf{P}\mathbf{X} = \mathbf{X}$ ("projection") and $\mathbf{M}\mathbf{X} = \mathbf{0}$ ("annihilator") <br>
    Also $\hat{\mathbf{y}} = \mathbf{P}\mathbf{y}$ ("hat") and $\mathbf{e} = \mathbf{M}\mathbf{y}$ ("residual maker") <br>

    <h5>5. Weak Law of Large Numbers, Central Limit Theorem</h5>

    Weak LLN: $\lim_{n\to\infty}\Pr\!\left(\,|\overline{X}_n-\mu| < \varepsilon\,\right) = 1$, where $\varepsilon > 0$<br>
    CLT: $\sqrt{n}\left(\bar{X}_n - \mu\right)\ \xrightarrow{d}\ \mathcal{N}\left(0,\sigma^2\right)$

    <h5>6. Slutksy's Theorem, Continuous Convergence Theorem, Asymptotic Equivalence Lemma</h5>

    Slutsky: If $X_n \xrightarrow{d} X$ and $Y_n \xrightarrow{p} c$ constant, then $X_n + Y_n \xrightarrow{d} X + c$ and $Y_nX_n \xrightarrow{d} cX$<br>
    CMT: If $X_n \xrightarrow{d} X$ and $g$ continuous, then $g(X_n) \xrightarrow{d} g(X)$ <br>
    Asymptotic equivalence: $X_n$ and $Y_n$ are asymptotically equivalent if $X_n - Y_n \xrightarrow{p} 0$
    
    <h5>7. Big Op, Little op, and the algebra of them</h5>

    Given $X_n$ is a RV, $a_n$ is constant, $\varepsilon >0$:<br>
    $X_n = O_p(a_n)$ means $P(|\frac{X_n}{a_n}| > \delta) < \epsilon, \forall n > N$<br>
    If $X_n = o_p(a_n)$, then $\frac{x_n}{a_n} = o_p(1)$, meaning $\lim_{n\to\infty} (P|\frac{X_n}{a_n}| \geq \epsilon) = 0, \forall\epsilon > 0$
    or in other words $\frac{X_n}{a_n} \xrightarrow{p} 0$<br>
    $o_p(a_n)$ implies $O_p(a_n)$ but not vice-versa.

    <h5>8. Delta Method</h5>

    Start with $X_n \sim N(\mu,\frac{\sigma^2}{n})$ and get $\frac{\sqrt{n}(X_n - \mu)}{\sigma} \sim N(0,1)$ <br>
    If $g$ is smooth then by Delta Method:<br>
    $\frac{\sqrt{n}(g(X_n) - g(\mu))}{|g'(\mu)|\sigma} \approx N(0,1)$ 
    and $g(X_n) \approx N\left(g(\mu), \frac{g'(\mu)^2\sigma^2}{n}\right)$<br>
    Note the approximation because $g(X_n)$ is an infinite sum.


    <h5>9. Frisch-Waugh Partialling Out</h5>

    Take $\mathbf{y}$ and two sets of regressors $\mathbf{X}_1$ and $\mathbf{X}_2$.
    FWL theorem claims that $\mathbf{b}_2$ will be equal after estimating in the following two forms:<br>
    $\mathbf{y} = \mathbf{X}_1 \mathbf{b}_1 + \mathbf{X}_2 \mathbf{b}_2 + \mathbf{e}$,
            then $\mathbf{b}_2 = (\mathbf{X}_2'\mathbf{X}_2)^{-1}\mathbf{X}_2'\mathbf{y}$ <br>
    $\mathbf{y}^* = \mathbf{X}^*_2 \mathbf{b}_2 + \mathbf{e}$,
        then $\mathbf{b}_2 = (\mathbf{X}^{* \prime}_2 \mathbf{X}^{*}_2 )^{-1} \mathbf{X}^{* \prime}_2 \mathbf{y}^*$, 
            where $\mathbf{X}^*_2 = \mathbf{M}_1 \mathbf{X}_2$ and $\mathbf{y}^* = \mathbf{M}_1 \mathbf{y}$<br>
    Explanation for notation in the second form:<br>
        i. Regress $\mathbf{y}$ on $\mathbf{X}_1$ and get the residuals $\mathbf{M}_1 \mathbf{y}$<br>
        ii. Regress $\mathbf{X}_2$ on $\mathbf{X}_1$ and get the residuals $\mathbf{M}_{1} \mathbf{X}_2$ <br>
        iii. Regress $\mathbf{M}_1 \mathbf{y}$ (residuals from i.) on $\mathbf{M}_1 \mathbf{X}_2$ (residuals from ii.), which will give $\mathbf{b}_2$ <br>


    <h5>10. For PD matrices A and B, A - B is PSD if and only if B^(-1) - A^(-1) is PSD: </h5>

    Sketch of the proof, "=>" direction:<br>
    Take arbitrary non-zero vector $x$, then $x' (A-B)x \geq 0$ and $A \geq B$<br>
    Then "divide" both sides by AB to get $A^{-1}AB^{-1} \geq A^{-1}BB^{-1}$ and then $B^{-1} \geq A^{-1}$ <br>
    Then $x' B^{-1} \geq x'A^{-1}$ and then $x' (B^{-1} - A^{-1})x \geq 0$<br>
    "<=" direction: replace $A$ with $B^{-1}$ and $B$ with $A^{-1}$ and repeat the same steps<br>

        """,
        unsafe_allow_html=True,
    )

    st.markdown(r"""<h4> Bonus (not Wooldridge)</h4>""", unsafe_allow_html=True)
    st.markdown(r"""11. **Concepts of Independence**""", unsafe_allow_html=True)

    with st.expander("Click to expand"):
        st.markdown(
            r"""
      
        Independence: $ P(X \cap Y) = P(X)P(Y) $ or $ P(X = x, Y = y) = P(X = x)P(Y = y) \; \forall x, y$. <br>
        Mean Independence: $ E(X | Y = y) = E(X) \; \forall y$.<br>
        Orthogonality: $\langle X, Y \rangle  = E[XY]=0$, but $\langle X, Y \rangle$ often defined as $Cov(X,Y)$, then
        $\langle X, Y \rangle = Cov(X,Y) = E[(X - E(X))(Y - E(Y))] = 0 $<br>
        Conditional Independence: $ P(X = x, Y = y | Z = z) = P(X = x | Z = z)P(Y = y | Z = z) \; \forall x, y, z$.<br>
        Linear Independence (in linear algebra):<br>
        Vectors $ \mathbf{v_1}, \mathbf{v_2}, ... \mathbf{v_n} $ are linearly independent if and only if no vector is zero and no vector can be expressed as a linear combination of the others.<br>
""",
            unsafe_allow_html=True,
        )

    st.markdown(
        r"""
        12. **Concepts of Convergence**<br>""",
        unsafe_allow_html=True,
    )

    with st.expander("Click to expand"):
        st.markdown(
            r"""
        Let $X$ be a random variable and $X_n$ a sequence of random variables.<br>
        i. Convergence in probability $X_n \xrightarrow{P} X$ if $ \lim_{n \to \infty} P(|X_n - X| > \epsilon) = 0 \; \forall \epsilon > 0$ <br>
        ii. Convergence almost surely: $X_n \xrightarrow{a.s.} X$ if $P(\lim_{n \to \infty} X_n = X) = 1$ <br>
        iii. Convergence in distribution: $X_n \xrightarrow{d} X$ if $\lim_{n \to \infty} F_{X_n}(x) = F_X(x) \; \forall x \in \mathbb{R}$ at which $F_X$ is continuous.<br>
        $ i. \implies ii. \implies iii.$ 

        """,
            unsafe_allow_html=True,
        )

pages/0_tools.py
justgri/phd_econometrics
import os
import random
import sys

import numpy as np
import pandas as pd
import src.scripts.plot_themes as thm
import src.scripts.utils as utl
import statsmodels.api as sm
import streamlit as st
from matplotlib import pyplot as plt
from scipy.stats import t
from st_pages import add_page_title

### PAGE CONFIGS ###

st.set_page_config(
    page_title="PhD Econometrics - OLS",
    page_icon="📈",
    layout="wide",
)
utl.local_css("src/styles/styles_pages.css")

random_seed = 0

### Data viz part
# Generate and display fixed data (data points, reg line, confidence interval)
# Show table with coefficients with their standard errors
# Then let user to choose underlying parameters for a, b, variance of errors and N
# Plot on the same chart and add coefficients to the table for comparison
# Bonus: show how coef SE depends on N and var(e) graphically

### Theory part
# show matrix notation, verify that sklearn gives same as multiplying by hand
# some experiments with projection matrices (P, M)?

# create one column with consistent width
s1, c01, s2 = utl.wide_col()

### PAGE START ###
# Dashboard header

with c01:
    st.title("Ordinary Least Squares Estimation")
    st.divider()
    st.header("1. Visualizing OLS estimates")

    st.markdown(
        r"""Suppose the true population relationship between $X$ and $y$ is defined by the slider values below.<br>
        You then draw a sample of size $n$ from that population and estimate OLS coefficients, $b_0$ and $b_1$.<br>
        Move the sliders to see how the sample estimates change with different population parameters.""",
        unsafe_allow_html=True,
    )
    st.latex(
        r"""
            y_i = \beta_0 + \beta_1x_i + \varepsilon_i \text{, where }  \varepsilon \sim N(0, \sigma^2) \\
            \hat{y_i} = b_0 + b_1 x_i
        """
    )


def gen_lin_data(b0, b1, sd, N, rseed):
    np.random.seed(rseed)
    # generate X
    K = 2
    X = np.random.uniform(-10, 10, (N, K - 1))
    X = np.column_stack((np.ones(N), X))

    # generate  error term
    eps = np.random.normal(0, sd, N)

    # y = xB + eps
    y = np.dot(X, np.array([b0, b1])) + eps

    # fit reg
    model = sm.OLS(y, X).fit()

    # get fitted values and CI
    predictions = model.get_prediction(X)
    y_hat = predictions.predicted_mean
    y_hat_se = predictions.se_mean

    # get 95% confidence interval
    ci = predictions.conf_int(alpha=0.05)  # 95% CI
    deg_freedom = X.shape[0] - X.shape[1]  # n - K

    # get error estimates
    e = y - y_hat
    s = np.sqrt(np.sum(e**2) / deg_freedom)

    return {
        "y": y,
        "x": X,
        "s": s,
        "y_hat": y_hat,
        "ci": ci,
        "model": model,
    }


slider_col, s1, chart_col = st.columns((0.8, 0.1, 1))

with slider_col:
    # Sliders
    b0_cust = st.slider(
        r"Intercept, $\beta_0$",
        min_value=-10.0,
        max_value=10.0,
        value=0.0,
        step=0.1,
    )
    b1_cust = st.slider(
        r"Slope, $\beta_1$", min_value=-5.0, max_value=5.0, value=0.0, step=0.1
    )
    var_cust = st.slider(
        r"Error SD, $\sqrt{var(\varepsilon)} = \sigma$",
        min_value=0.1,
        max_value=20.0,
        value=10.0,
        step=0.1,
    )

    n_cust = st.slider(
        "Sample size, $n$",
        min_value=10,
        max_value=1000,
        value=500,
        step=10,
    )

custom_data = gen_lin_data(b0_cust, b1_cust, var_cust, n_cust, random_seed)


def plot_ols(data_custom, b0, b1):
    fig, ax = plt.subplots(figsize=(4, 4))
    plt.subplots_adjust(left=0)  # remove margin

    # Sample data
    ax.scatter(
        data_custom["x"][:, 1],
        data_custom["y"],
        # label="Custom Data",
        color=thm.cols_set1_plt[1],
        alpha=0.5,
    )

    # Fitted line
    b0_s = data_custom["model"].params[0]
    b1_s = data_custom["model"].params[1]

    if b1_s >= 0:
        label_s = rf"$\hat{{y}} = {b0_s:.2f} + {b1_s:.2f}x$"
    else:
        label_s = rf"$\hat{{y}} = {b0_s:.2f} - {-b1_s:.2f}x$"

    ax.plot(
        data_custom["x"][:, 1],
        data_custom["y_hat"],
        label=label_s,
        color=thm.cols_set1_plt[4],
        linewidth=2,
    )

    sorted_indices = np.argsort(data_custom["x"][:, 1])
    sorted_x = data_custom["x"][:, 1][sorted_indices]
    sorted_ci_lower = data_custom["ci"][:, 0][sorted_indices]
    sorted_ci_upper = data_custom["ci"][:, 1][sorted_indices]

    ax.fill_between(
        sorted_x,
        sorted_ci_lower,
        sorted_ci_upper,
        color=thm.cols_set1_plt[4],
        alpha=0.3,
        label="95% Confidence Interval",
    )

    plt.xlim([-11, 11])
    plt.ylim([-50, 50])
    plt.xlabel("X", fontweight="bold")
    plt.ylabel("Y", fontweight="bold")
    ax.yaxis.set_label_coords(-0.08, 0.5)
    plt.legend(loc="upper left", fontsize="small")

    return fig


def create_summary(data):
    coefficients = pd.DataFrame(
        {
            "Coefficient": [r"Intercept", "Slope"],
            "Population Parameters": [b0_cust, b1_cust],
            "Sample Estimates": [
                data["model"].params[0],
                data["model"].params[1],
            ],
            "Standard Errors": [
                data["model"].bse[0],
                data["model"].bse[1],
            ],
        }
    )

    # Apply formatting to the "True Pop" and "Estimate" columns
    coefficients[
        ["Population Parameters", "Sample Estimates", "Standard Errors"]
    ] = coefficients[
        ["Population Parameters", "Sample Estimates", "Standard Errors"]
    ].applymap(
        lambda x: f"{x:.2f}"
    )

    return coefficients


with slider_col:
    if st.button("Resample data", type="primary"):
        random_seed = random.randint(0, 10000)
        custom_data = gen_lin_data(b0_cust, b1_cust, var_cust, n_cust, random_seed)


coefficients = create_summary(custom_data)

with chart_col:
    chart_col.pyplot(plot_ols(custom_data, b0_cust, b1_cust), use_container_width=True)


# CSS styles for the table (center and header)
table_styler = [
    {
        "selector": "th",  # Apply to header cells
        "props": [("background-color", "lightblue")],
    },
    {
        "selector": "td",  # Apply to data cells
        "props": [
            ("font-size", "20px"),
            ("text-align", "center"),
            ("background-color", "white"),
        ],
    },
]

# Apply CSS styles to the DataFrame
styled_coefficients = coefficients.style.set_table_styles(table_styler)

# Create a centered and styled HTML representation of the DataFrame
styled_table = styled_coefficients.hide(axis="index").to_html()

# Define custom CSS to style the table and center it
table_css = """
<style>
    table {
        margin: 0 auto;
        text-align: center;
    }
</style>
"""

s0, c02, s1 = utl.narrow_col()

with c02:
    st.write("")
    # display table and plot
    st.write(
        f"{table_css}{styled_table}",
        unsafe_allow_html=True,
    )
    st.write("")

    st.latex(
        f"n= {n_cust} ,"
        + f"R^2 = {custom_data['model'].rsquared:.2f}"
        + r", \hat{\sigma}"
        + f"= {custom_data['s']:.2f}"
    )

s0, c03, s1 = utl.wide_col()

with c03:
    st.markdown("### Interesting takeaways")

    with st.expander("Click to expand."):
        st.markdown(
            r"""                
        1. Variance of $\mathbf{b}$ (and thus the standard error) does not depend on population $\beta$. <br>
        It depends on variance of errors $s^2$ (and thus $\sigma^2)$, $N-k$, and $X'X$.<br>
        Note that higher variance of $X$ leads to a lower variance of $\mathbf{b}$, which is intuitive because you cover a wider range of $X$s.<br>
        
            $\widehat{var(\mathbf{b}| \mathbf{X})} \equiv s^2 (X'X)^{-1} = \frac{\mathbf{e'e}}{n - K} (X'X)^{-1}$
    """,
            unsafe_allow_html=True,
        )

        # ADD FORMULAS FOR CONFIDENCE INTERVALS
        # check against statsmodels link below
        # Greene Ch 4.8 (p. 86)
        # https://www.statsmodels.org/0.9.0/_modules/statsmodels/regression/_prediction.html#PredictionResults
    #     st.markdown(
    #         r"""
    #     e. Confidence interval is visually more sensitive to sample size than to error variance, i.e.,
    #     if $n$ is large enough, even for high $\sigma^2$, the confidence interval is small because it depends only on $s^2$ and not on $\sigma^2$:<br>
    #     Also, CI is wider when further away from mean.<br>
    # """,
    #         unsafe_allow_html=True,
    #     )


with c03:
    st.header("2. OLS assumptions")
    st.write("Greene Chapter 4, Table 4.1:")
    st.markdown(
        r"""
        **A1. Linearity:** $\mathbf{y = X \beta + \varepsilon}$ <br>
        **A2. Full rank:** $\mathbf{X}$ is an $ n \times K$ matrix with rank $K$ (full column rank) <br>
        **A3. Exogeneity of the independent variables:** $E[\varepsilon_i | \mathbf{X}] = 0$ <br>
        **A4. Homoscedasticity and nonautocrrelation:** $E[\mathbf{\varepsilon \varepsilon'} | \mathbf{X}] = \sigma^2 \mathbf{I}_n$ <br>
        **A5. Stochastic or nonstochastic data:** $\mathbf{X}$ may be fixed or random.<br>
        **A6. Normal distribution:** $\varepsilon | \mathbf{X} \sim N(0, \sigma^2 \mathbf{I}_n)$ <br>
        """,
        unsafe_allow_html=True,
    )

s0, c04, s1 = utl.wide_col()

with c04:
    st.header("3. Theory with code")

    def tabs_code_theory():
        return st.tabs(["Theory", "Code numpy", "Code statsmodels"])

    st.markdown(
        "#### Solving for OLS coefficients",
        unsafe_allow_html=True,
    )
    # OLS Theory vs Code
    f1_c1, f1_c2, f1_c3 = tabs_code_theory()
    with f1_c1:
        st.markdown(
            r"""
            Let $\mathbf{x}_i$ and $\beta$ be $K \times 1$ vectors, and $\mathbf{X}$ be $n \times K$ matrix of $\mathbf{x}_i$<br>
            $y_i = \mathbf{x}_i' \beta + \varepsilon_i, (i = 1, 2, ..., n)$<br>
            $\mathbf{y = X \beta + \varepsilon}$<br>
            <br>
            Let $\mathbf{e \equiv y - Xb}$ and solve for $\mathbf{b}$ that minimizes $\mathbf{e'e}$<br>
            Solution to this minimization problem is:<br>
            $\mathbf{b = (X'X)^{-1}X'y}$, which we call the OLS estimate for $\beta$<br>
            And then we can get predicted ("fitted") $y$ values:
            $\mathbf{\hat{y} = X b}$ <br>
            <br>
            $\mathbf{P_{n \times n} \equiv X(X'X)^{-1}X'}$ is the "projection" or "hat" matrix<br>
            $\mathbf{PX = P}$ and $\mathbf{Py = \hat{y}}$<br>
            <br>
            $\mathbf{M_{n \times n} \equiv I_n - P}$ is the "annihilator" or "residual maker" matrix<br>
            $\mathbf{MX = 0}$ and $\mathbf{My = e}$<br>
            

    """,
            unsafe_allow_html=True,
        )

    with f1_c2:
        ols_code = """
        import numpy as np
        # Define initial parameters
        n = 1000 # sample size
        K = 2 # number of variables and coefficients
        sigma = 1 # standard deviation of error term
        X = np.random.uniform(-10, 10, (n, K - 1)) # generate random X
        X = np.column_stack((np.ones(n), X)) # add constant
        beta = np.random.uniform(-5, 5, K) # generate random true beta
        eps = np.random.normal(0, sigma, n) # generate random errors (disturbances)

        # Generate true y
        y = np.dot(X, beta) + eps

        # Solve for OLS coefficients
        b_ols = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(y)
        y_hat = X.dot(b_ols)
        e = y - y_hat

        # Get projection and residual maker matrices
        P = X.dot(np.linalg.inv(X.T.dot(X))).dot(X.T)
        M = np.eye(n) - P

        """
        st.code(ols_code, language="python")

    with f1_c3:
        ols_code_b = """
        import statsmodels.api as sm
        # Generate data with numpy same as before
        # ...

        # Fit OLS model with statsmodels
        model = sm.OLS(y, X).fit()
        b_ols = model.params
        y_hat = model.predict(X)
        e = model.resid

        # statsmodels does not have built-in P and M matrices
        # since you can pull necessary estimates from model.get_prediction(X)

        """

        st.code(ols_code_b, language="python")

    st.divider()

    ### OLS Additional Concepts
    st.markdown(
        "#### OLS Errors",
        unsafe_allow_html=True,
    )

    f2_c1, f2_c2, f2_c3 = tabs_code_theory()
    with f2_c1:
        # Prediction error (MSPE) - Greene p. 86
        # review the following sources
        # https://dept.stat.lsa.umich.edu/~kshedden/Courses/Regression_Notes/prediction.pdf
        # http://web.vu.lt/mif/a.buteikis/wp-content/uploads/PE_Book/3-7-UnivarPredict.html

        st.markdown(
            r"""
            OLS estimate for $\sigma^2$ aka Standard Error of the Regression (SER):<br>
            $\hat{\sigma}^2 = s^2 \equiv \frac{SSE}{n-K} = \frac{\mathbf{e'e}}{n-K}$<br>
            $\hat{\sigma} = \sqrt{s^2} = s$<br>
        
            Conditional variance of $\mathbf{b}$ and standard error of $b_k$:<br>
            $Var(\mathbf{b|X}) = s^2(X'X)^{-1}$<br>
            $SE(b_k) = \{[s^2(X'X)^{-1}]_{kk}\}^{1/2}$<br>
            (square root of *k*th diagonal element of the variance matrix)

            Prediction error and its variance - take a new vector $\mathbf{x}^0$ for which you want to predict dependent variable $\mathbf{y}^0$:<br>
            $\mathbf{y}^0 = \mathbf{x}^{0\prime} \mathbf{b} + \varepsilon^0 = \mathbf{x}^{0\prime} \mathbf{b}$ (by A3)<br>
            $e^0 = \hat{y}^0 - y^0 = \mathbf{(b - \beta)'x}^0 - \varepsilon^0$<br>
            $Var[e^0 | \mathbf{X, x}^0] = \sigma^2 + \mathbf{x}^{0\prime} [\sigma^2 (X'X)^{-1}]\mathbf{x}^0$<br>
            
            Prediction variance can be estimated by using $s^2$ instead of $\sigma^2$, thus the standard error of prediction error is:<br>
            $se(e^0) = \sqrt{s^2 + \mathbf{x}^{0\prime} [s^2 (X'X)^{-1}]\mathbf{x}^0}$<br>
            
            Prediction interval = $\hat{y}^0 \pm t_{\alpha/2, n-K} se(e^0)$<br>

            Sampling error:<br>
            $\mathbf{b} - \beta = \mathbf{(X'X)^{-1}X'y - \beta} = \mathbf{(X'X)^{-1}X' \varepsilon}$ (sampling error)<br>
            (only relevant if $\beta$ is known)

         """,
            unsafe_allow_html=True,
        )

    with f2_c2:
        ols_code_err = """
        import numpy as np

        # OLS estimate for sigma^2 and sigma
        s_sq = e.dot(e) / (n - K)
        s = np.sqrt(s_sq)

        # Sampling error (only b/c we know true beta values)
        sampling_error = b_ols - beta

        # Conditional variance and standard errors of beta
        var_b = s_sq * np.linalg.inv(X.T.dot(X))
        SE_b = np.sqrt(np.diag(var_b))
        """
        st.code(ols_code_err, language="python")

        # Prediction confidence intervals - double check
        # deg_freedom = X.shape[0] - X.shape[1]  # n - K
        # alpha = 0.05
        # t_score = t.ppf(1 - alpha/2, deg_freedom)
        # ci = np.column_stack(
        #     (y_hat - t_score * y_hat_se, y_hat + t_score * y_hat_se)
        # )

    with f2_c3:
        ols_code_err_b = """
        import statsmodels.api as sm
        model = sm.OLS(y, X).fit()
        
        # OLS estimate for sigma^2 and sigma
        s_sq = model.mse_resid
        s = model.mse_resid ** 0.5
        
        # Sampling error (only b/c we know true beta values)
        sampling_error = model.params - beta
        
        # Conditional variance and standard errors of beta
        var_b = model.cov_params()
        SE_b = model.bse

        """
        # Prediction confidence intervals - need to add
        # predictions = model.get_prediction(X)
        # y_hat_se = predictions.se_mean
        # ci = predictions.conf_int(alpha=0.05)  # 95% CI
        # st, data, ss2 = summary_table(re, alpha=0.05)
        # fittedvalues = data[:, 2]
        # predict_mean_se = data[:, 3]
        # predict_mean_ci_low, predict_mean_ci_upp = data[:, 4:6].T # which is same as ci
        # predict_ci_low, predict_ci_upp = data[:, 6:8].T

        st.code(ols_code_err_b, language="python")


s0, c05, s1 = utl.wide_col()
with c05:
    st.header("4. More details")
    st.write("Check out Matteo Courthoud's website:")
    st.link_button(
        "OLS Algebra",
        "https://matteocourthoud.github.io/course/metrics/05_ols_algebra/",
        type="primary",
    )

    st.header("5. Proofs to remember")
    sst_proof = "https://stats.stackexchange.com/questions/207841/why-is-sst-sse-ssr-one-variable-linear-regression/401299#401299"

    with st.expander("OLS estimate is BLUE"):
        st.link_button(
            "See Matteo Courthoud's website",
            "https://matteocourthoud.github.io/course/metrics/05_ols_algebra/#blue",
            type="secondary",
        )

    with st.expander("Transformed variables"):
        st.markdown(
            r"""
            Proof from Greene THEOREM 3.8:<br>
            Let $\mathbf{b}$ be the OLS coefficients from regressing $\mathbf{y}$ on $\mathbf{X}$.<br>
            Regress $\mathbf{y}$ on transformed variables $\mathbf{Z = XP}$, where $\mathbf{P}$ transforms columns of $\mathbf{X}$.<br>
            Then OLS coefficients are:<br>
            $\mathbf{d = (Z'Z)^{-1}Z'y = [(XP)'(XP)]^{-1}(XP)'y}$ $\mathbf{= [P'X'XP]^{-1}P'X'y =}$
            $\mathbf{P^{-1}(X'X)^{-1}P'^{-1}P'X'y}$ $\mathbf{= P^{-1}(X'X)^{-1}X'y = P^{-1}b}$<br>
            The vector of residuals is:<br>
            $\mathbf{u = y - Z(P^{-1}b) = y - XP(P^{-1}b) = y - Xb = e}$<br>
            Since residuals are identical and $\mathbf{y'M^0y}$ is unchanged, $R^2$ is also identical in two regressions.
            """,
            unsafe_allow_html=True,
        )

pages/1_ols.py
justgri/phd_econometrics
import os
import random
import sys

import numpy as np
import pandas as pd
import src.scripts.plot_themes as thm
import src.scripts.utils as utl
import statsmodels.api as sm
import streamlit as st
from matplotlib import pyplot as plt
from scipy.stats import t
from st_pages import add_page_title

### PAGE CONFIGS ###

st.set_page_config(
    page_title="PhD Econometrics - Fit",
    page_icon="📈",
    layout="wide",
)

utl.local_css("src/styles/styles_pages.css")

random_seed = 0

## Data viz part
# How R-sq is calculated
# How R-sq and other measures depend on sample size, error variance, and true parameters
# Also how much R-sq, AIC, and BIC are penalized by adding more regressors?

## Theory part
# TBD

# create one column with consistent width
s1, c01, s2 = utl.wide_col()

### PAGE START ###
# Dashboard header

with c01:
    st.title("Model fit measures")
    st.divider()
    st.header("1. Visualizing R-squared")

    st.markdown(
        r"""Suppose the true population relationship between $X$ and $y$ is linear and defined by the slider values below.
        Move the sliders to see how the model fit measures change with different population parameters.""",
        unsafe_allow_html=True,
    )


def gen_lin_data(b0, b1, sd, N, rseed):
    np.random.seed(rseed)
    # fix X values
    X_fix = [-2, -0.5, 2, 0.5, -1.25, 1.25, 0]
    # X_pos_ref = [0, 2, 6, 4, 1, 5, 3]

    K = 2
    X = np.array(X_fix[:N])
    X = np.column_stack((np.ones(len(X)), X))

    # generate  error term
    eps_fix = np.array([-1 * sd, sd, -sd, -sd, -1 * sd, 1 / 2 * sd, -sd])
    eps = np.array(eps_fix[:N])
    # eps = np.random.normal(0, sd, N)

    # y = xB + eps
    y = np.dot(X, np.array([b0, b1])) + eps

    # fit reg
    model = sm.OLS(y, X).fit()

    # get fitted values and CI
    predictions = model.get_prediction(X)
    y_hat = predictions.predicted_mean
    y_hat_se = predictions.se_mean

    # get error parameters
    e = y - y_hat
    deg_freedom = X.shape[0] - X.shape[1]  # n - K
    s = np.sqrt(np.sum(e**2) / deg_freedom)

    # calculate R^2 manually
    y_bar = np.mean(y)
    ss_tot = np.sum((y - y_bar) ** 2)
    ss_res = np.sum((y - y_hat) ** 2)
    r_squared = 1 - ss_res / ss_tot

    return {
        "y": y,
        "x": X,
        "s": s,
        "y_hat": y_hat,
        "model": model,
    }


slider_col, s1, chart_col = st.columns((0.8, 0.1, 1))

with slider_col:
    # Sliders
    b0_cust = st.slider(
        r"Intercept, $\beta_0$",
        min_value=-1.5,
        max_value=1.5,
        value=0.0,
        step=0.1,
    )
    b1_cust = st.slider(
        r"Slope, $\beta_1$", min_value=-1.5, max_value=1.5, value=1.0, step=0.1
    )
    var_cust = st.slider(
        r"Error SD, $\sqrt{var(\varepsilon)} = \sigma$",
        min_value=0.1,
        max_value=3.0,
        value=2.0,
        step=0.1,
    )

    n_cust = st.slider(
        "Points displayed, $n$",
        min_value=3,
        max_value=7,
        value=4,
        step=1,
    )

custom_data = gen_lin_data(b0_cust, b1_cust, var_cust, n_cust, random_seed)


def plot_ols(data_custom, b0, b1):
    fig, ax = plt.subplots(figsize=(4, 4))
    plt.subplots_adjust(left=0)  # remove margin
    # Define the upper and lower bounds for axes
    y_ub = 6
    y_lb = -y_ub
    x_ub = 3
    x_lb = -x_ub

    ### SAMPLE DATA
    ax.scatter(
        data_custom["x"][:, 1],
        data_custom["y"],
        # label="Custom Data",
        color=thm.cols_set1_plt[1],
        alpha=0.9,
        zorder=10,
        edgecolors="none",
        s=60,
    )

    ### REG LINE
    # Regression line
    b0_s = data_custom["model"].params[0]
    b1_s = data_custom["model"].params[1]

    if b1_s >= 0:
        label_s = rf"$\hat{{y}}_i = {b0_s:.2f} + {b1_s:.2f}x_i$"
    else:
        label_s = rf"$\hat{{y}}_i = {b0_s:.2f} - {-b1_s:.2f}x_i$"

    ax.plot(
        data_custom["x"][:, 1],
        data_custom["y_hat"],
        label=label_s,
        color=thm.cols_set1_plt[4],
        zorder=9,
        linewidth=2.5,
    )

    ### Y BAR
    y_mean = np.mean(data_custom["y"])
    plt.axhline(y=y_mean, color=thm.cols_set1_plt[2], linestyle="-", linewidth=2.5)

    # Y BAR text
    y_offset = 0.35 if b1_s > 0 else -0.45
    ax.text(
        x_lb + 0.1,
        y_mean + y_offset,
        r"$\bar{y}=$" + f"{y_mean:.2f}",
        color="green",
        verticalalignment="center",
    )

    ### REG RESIDUALS (for residual sum of squares)
    res_offset = 0.05
    has_lab_res = False
    for xi, yi, yhati in zip(
        data_custom["x"][:, 1], data_custom["y"], data_custom["y_hat"]
    ):
        if not has_lab_res:
            lab_res = r"$y_i - \hat{y_i}$ for SSE"
            has_lab_res = True
        else:
            lab_res = None  # No label for other lines

        plt.plot(
            [xi - res_offset, xi - res_offset],
            [yi, yhati],
            color=thm.cols_set1_plt[4],
            linestyle="--",
            alpha=0.9,
            label=lab_res,
            linewidth=2,
        )

    ### DEVIATIONS FROM THE MEAN (for total sum of squares)
    has_lab_mean = False
    # Add horizontal lines from each data point to the mean of Y (Total Sum of Squares)
    for xi, yi in zip(data_custom["x"][:, 1], data_custom["y"]):
        if not has_lab_mean:
            lab_mean = r"$y_i - \bar{y}\;$ for TSS"
            has_lab_mean = True
        else:
            lab_mean = None  # No label for other lines

        plt.plot(
            [xi + res_offset, xi + res_offset],
            [yi, y_mean],
            color=thm.cols_set1_plt[2],
            linestyle=":",
            alpha=0.9,
            label=lab_mean,
            linewidth=2,
        )

    # Add R-sq label
    r_squared = data_custom["model"].rsquared

    # Add R^2 to the top right corner
    rsq_offset = 0.1
    ax.text(
        x_ub - rsq_offset,
        y_ub - rsq_offset - 0.1,
        f"$R^2 = {r_squared:.2f}$",
        horizontalalignment="right",
        verticalalignment="top",
        fontsize=10,
    )

    plt.xlim([x_lb, x_ub])
    plt.ylim([y_lb, y_ub])
    plt.xlabel("X", fontweight="bold")
    plt.ylabel("Y", fontweight="bold")
    ax.yaxis.set_label_coords(-0.08, 0.5)

    legend_loc = "upper left" if data_custom["y_hat"][0] < 3.25 else "lower left"
    plt.legend(loc=legend_loc, fontsize="small")

    return fig


with chart_col:
    chart_col.pyplot(plot_ols(custom_data, b0_cust, b1_cust), use_container_width=True)


def generate_html_table(model):
    # Extract statistics
    ess = np.sum(model.resid**2)
    tss = np.sum((model.model.endog - np.mean(model.model.endog)) ** 2)
    r2 = model.rsquared
    r2_adj = model.rsquared_adj

    likelihood = model.llf
    aic = model.aic
    bic = model.bic
    N = model.nobs

    data = [
        ("SSE (Error Sum of Squares)", ess),
        ("TSS (Total Sum of Squares)", tss),
        ("R-squared", r2),
        ("Adjusted R-squared", r2_adj),
        # ("Likelihood", likelihood),
        # ("AIC", aic),
        # ("BIC", bic),
    ]

    html_string = """<table class="table" border="1">
                    <thead>
                        <tr>
                            <th>Name</th>
                            <th>Value</th>
                        </tr>
                    </thead>
                    <tbody>"""

    for name, value in data:
        html_string += f"<tr><td>{name}</td><td>{value:.2f}</td></tr>"

    html_string += "</tbody></table>"

    return html_string


table_html_string = generate_html_table(custom_data["model"])

_, formula_col, s1, table_col, _ = st.columns((0.2, 1, 0.1, 1, 0.2))

with table_col:
    st.markdown(table_html_string, unsafe_allow_html=True)

with formula_col:
    st.markdown(r"""<h5>Calculations</h5>""", unsafe_allow_html=True)
    st.markdown(
        r"""
        $SSE = \sum_{i=1}^n (y_i-\hat{y_i})^2 = \mathbf{e'e}$<br>
        $TSS = \sum_{i=1}^n (y_i-\bar{y_i})^2 = \mathbf{y'M^0y}$<br>
        $R^2 = 1 - \frac{SSE}{SST} = 1- \mathbf{\frac{e'e}{y'M^0y}}$<br>
        $\bar{R}^2 = 1 - \frac{n - 1}{n - K} (1 - R^2)$<br>
""",
        unsafe_allow_html=True,
    )
_, nb_col, _ = st.columns((0.1, 1, 0.1))
nb_col.markdown(
    r"NB: Adjusted R-squared $\bar{R}^2$ is heavily penalized here, because the sample size is very small.",
    unsafe_allow_html=True,
)


s0, c03, s1 = utl.wide_col()

with c03:
    st.markdown("### Interesting takeaways")

    with st.expander("Click to expand."):
        st.markdown(
            r"""
        1. $R^2$ is expected to be 0 if $\beta_1=0$ or if $X_i = \bar{X}$.<br>

            $R^2 = \frac{ (\hat{y} - \bar{y})' (\hat{y} - \bar{y}) }{ (y - \bar{y})' (y - \bar{y}) }
            = \frac{\sum_{i=1}^{N} (\hat{y}_i - \bar{y})^2}{\sum_{i=1}^{N} (y_i - \bar{y})^2}
            = \frac{\sum_{i=1}^{N} (\hat{\beta_1} (X_i - \bar{X}))^2}{\sum_{i=1}^{N} (y_i - \bar{y})^2}$
            , because $\hat{\beta_0} = \bar{y} - \hat{\beta_1}\bar{X}$

        2. $R^2$ is independent of the size of the intercept $b_0$ but is dependent on the size of the slope $b_1$.
        
            """,
            unsafe_allow_html=True,
        )

        # Explain why R^2 increases with \beta - what's the interpretation?

        # ADD FORMULAS FOR CONFIDENCE INTERVALS

        # related - show that R and CI are related?
        # check against statsmodels link below
        # https://www.statsmodels.org/0.9.0/_modules/statsmodels/regression/_prediction.html#PredictionResults
    #     st.markdown(
    #         r"""
    #     e. Confidence interval is visually more sensitive to sample size than to error variance, i.e.,
    #     if $n$ is large enough, even for high $\sigma^2$, the confidence interval is small because it depends only on $s^2$ and not on $\sigma^2$:<br>
    #     Also, CI is wider when further away from mean.<br>
    # """,
    #         unsafe_allow_html=True,
    #     )


with c03:
    st.header("2. Interpretation")
    st.markdown(r"""<h5>R-squared indicates:</h5>""", unsafe_allow_html=True)
    st.markdown(
        r"""
    1. "Whether variation in $x$ is a good predictor of variation in $y$" - useful for prediction accuracy (G Ch 3.5)<br>
    2. Proportion of total variation in $y$ that is accounted by variation in $x$ - useful for model comparison<br>
        NB: The dependent variable has to be of the same transformation (e.g., linear or log) for comparison.       
    3. How close the data are to the estimated model - useful to get a sense of over/underfitting
""",
        unsafe_allow_html=True,
    )

    st.markdown(r"""<h5>R-squared does NOT indicate:</h5>""", unsafe_allow_html=True)
    st.markdown(
        r"""
        1. Causality
        2. Unbiasedness of the coefficients
        3. Appropriateness of the model
        4. How reasonable OOS predictions are (e.g., high R-sq doesn't prevent forecasting of negative wages)

        NB: In general, **interpretation of the R-squared values is highly dependent on the context**. 
        """,
        unsafe_allow_html=True,
    )

    st.markdown(
        r"""<h5>R-sq and R-sq adj. vs AIC and BIC</h5>""",
        unsafe_allow_html=True,
    )

    st.markdown(
        r"""Since R-sq increases by including more regressors, we need a measure that penalizes for adding regressors.
        It has been argued that the adjusted R-sq doesn't penalize heavily enough, thus AIC and BIC measures have been proposed.(Greene p. 47)<br>""",
        unsafe_allow_html=True,
    )
    st.markdown(
        r"""
        "Choosing a model based on the lowest AIC is logically the same as using $\bar{R}^2$ in the linear model, nonstatistical, albeit widely accepted.
        The AIC and BIC are information criteria, not fit measures as such." (Greene, p. 561)<br>
        APC has a direct relationship to $R^2$.<br>
        """,
        unsafe_allow_html=True,
    )

s0, c04, s1 = utl.wide_col()

with c04:
    st.header("3. Theory with code")

    def tabs_code_theory():
        return st.tabs(["Theory", "Code statsmodels", "Code numpy"])

    ### Error sums
    st.markdown(
        "#### Error sums",
        unsafe_allow_html=True,
    )
    st.markdown(
        """NB: Hayashi and Greene classically disagree on notation for the sum of squared residuals (SSE or SSR), so I'll follow Greene.""",
        unsafe_allow_html=True,
    )

    f2_c1, f2_c2, f2_c3 = tabs_code_theory()
    with f2_c1:
        st.markdown(
            r"""
            Error Sum of Squares (SSE) aka Sum of Squared Residuals (SSR or RSS, hence confusion):<br>
            $SSE = \sum_{i=1}^n (y_i-\hat{y_i})^2 = \sum_{i=1}^n (e_i)^2 =  \mathbf{e'e = \varepsilon' M \varepsilon}$<br>
            (this is SSR according to Hayashi)<br>

            Regression sum of squares (SSR) aka Explained Sum of Squares (ESS):<br>
            $SSR = \sum_{i=1}^n (\hat{y_i} - \bar{y})^2 = \sum_{i=1}^n (\hat{y_i} - \bar{\hat{y}})^2$<br>
            $SSR =  \mathbf{b'X'M^0Xb}$, where $\mathbf{M^0}$ is the centering matrix<br>

            Total sum of squares (SST) aka Total Variation:<br>
            $SST = \sum_{i=1}^n (y_i-\bar{y_i})^2 = \sum_{i=1}^n (\hat{y_i} - \bar{y})^2 + \sum_{i=1}^n (e_i)^2$ <br>
            $SST = \mathbf{y'M^0y = b'X'M^0Xb + e'e = SSR + SSE}$<br>
         """,
            unsafe_allow_html=True,
        )

    with f2_c2:
        ols_code_err_b = """
        import statsmodels.api as sm
        model = sm.OLS(y, X).fit()

        # Sum of squared errors
        SSE = model.ssr # this is SSE according to Greene

        # Regression sum of squares
        SSR = model.ess
        
        # Total sum of squares
        SST = SSE + SSR
        """

        st.code(ols_code_err_b, language="python")

    with f2_c3:
        ols_code_err = """
        import numpy as np

        # Sum of squared errors
        SSE = e.dot(e)
        
        # Regression sum of squares
        y_hat_centered = y_hat - np.mean(y_hat)
        SSR = y_hat_centered.dot(y_hat_centered)

        # Total sum of squares
        y_centered = y - np.mean(y)
        SST = y_centered.dot(y_centered)
        """
        st.code(ols_code_err, language="python")

    st.divider()

    st.markdown("#### Model fit and selection measures")
    st.markdown(
        r"""
        NB: $R^2$ definition below requires a constant term to be included in the model.<br>
        """,
        unsafe_allow_html=True,
    )

    f3_c1, f3_c2, f3_c3 = tabs_code_theory()

    # Sources for AIC and BIC
    sas_source = "https://documentation.sas.com/doc/en/vfcdc/8.5/vfug/p0uawamu7dmtc2n1cllfwajyvlko.htm"
    stata_source = "https://www.stata.com/manuals13/restatic.pdf"
    stack_ex = "https://stats.stackexchange.com/questions/490056/aic-bic-formula-wrong-in-james-witten"

    with f3_c1:
        st.markdown(
            r"""          
            R-sq, Adjusted R-sq, and Pseudo R-sq:<br>
            $R^2 = \frac{SSR}{SST} = \frac{SST - SSE}{SST} = 1 - \frac{SSE}{SST}= 1- \mathbf{\frac{e'e}{y'M^0y}}$<br>
            $\bar{R}^2 = 1 - \frac{n - 1}{n - K} (1 - R^2)$<br>
            McFadden Pseudo  $R^2 = 1 - \frac{\text{ln} L}{\text{ln} L_0} = \frac{-\text{ln}(1-R^2)}{1+\text{ln}(2\pi) + \text{ln}(s_y^2)}$<br>
            
            Amemiya's Prediction Criterion (APC):<br>
            $APC=\frac{SSE}{n-K}(1+\frac{K}{n}) = SSE \frac{n+K}{n-K}$<br>

            AIC and BIC for OLS, when error variance is known (Greene p. 47):<br>
            $AIC = \text{ln}(\frac{SSE}{n}) + \frac{2K}{n}$<br>
            $BIC = \text{ln}(\frac{SSE}{n}) + \frac{\text{ln}(n) K}{n}$<br>
            
            AIC and BIC are more often calculated for any MLE as follows (Greene p. 561):<br>
            $AIC = -2 \text{ln}(L)+2K$<br>
            $BIC = -2 \text{ln}(L) + \text{ln}(n) K  $<br>
            
            In OLS, SSE is proportional to log-likelihood, so the two formulas would lead to the same model selection.<br>
            NB: Even for OLS, Python *statsmodels*, STATA *estat ic*, and R *lm* use the latter definition, whereas SAS uses the former multiplied by $n$.
            """,
            unsafe_allow_html=True,
        )

    with f3_c2:
        r2_code_built_in = """
        import statsmodels.api as sm
        import numpy as np

        model = sm.OLS(y, X).fit()

        # R-sq and R-sq adjusted
        r_sq = model.rsquared
        r_sq_adj = model.rsquared_adj

        # Pseudo R-sq
        ln_L = model.llf
        model_constant = sm.OLS(y, np.ones(n)).fit()
        ln_L_0 = model_constant.llf
        pseudo_r_sq = 1 - ln_L / ln_L_0

        # Amemiya's Prediction Criterion - no built-in module
        APC = (e.dot(e) / n) * (n + K) / (n - K)

        # AIC and BIC
        AIC = model.aic
        BIC = model.bic

        """

        st.code(r2_code_built_in, language="python")

    with f3_c3:
        r2_code = """
        import numpy as np
        # R-sq and R-sq adjusted
        y_centered = y - np.mean(y)
        r_sq = 1 - e.dot(e) / y_centered.dot(y_centered)
        r_sq_adj = 1 - ((n - 1) / (n - K)) * (1 - r_sq)

        # Pseudo R-sq
        var_y = np.var(y)
        pseudo_r_sq = (-1 * np.log(1 - r_sq) / (1 + np.log(2 * np.pi) + np.log(var_y)))
                
        # Amemiya's Prediction Criterion
        APC = (e.dot(e) / n) * (n + K) / (n - K)

        # AIC and BIC, first get log likelihood
        ln_L = (-n / 2) * (1 + np.log(2 * np.pi) + np.log(SSE / n))
        AIC = -2 * ln_L + 2 * K
        BIC = -2 * ln_L + K * np.log(n)
        """
        st.code(r2_code, language="python")

s0, c05, s1 = utl.wide_col()
with c05:
    st.header("4. More details")
    st.write("Check out Andrius Buteikis' book:")
    st.link_button(
        "Goodness-of-Fit",
        "http://web.vu.lt/mif/a.buteikis/wp-content/uploads/PE_Book/3-8-UnivarGoF.html",
        type="primary",
    )

    st.header("5. Proofs to remember")
    sst_proof = "https://stats.stackexchange.com/questions/207841/why-is-sst-sse-ssr-one-variable-linear-regression/401299#401299"

    with st.expander("SST = SSR + SSE"):
        st.markdown(
            rf"Proof from Greene Section 3.5 (also see [Stack Exchange]({sst_proof})):<br>"
            + r"""
                $y_i - \bar{y} = \mathbf{x}_i'\mathbf{b} + e_i$<br>
                $y_i - \bar{y} = \hat{y}_i - \bar{y} + e_i = (\mathbf{x}_i - \mathbf{\bar{x}})'\mathbf{b} + e_i$<br>
                $\mathbf{M^0y= M^0Xb + M^0e}$<br>
                $SST = \mathbf{y'M^0y = b'X'M^0Xb + e'e} = SSR + SSE$<br>
                (need to expand between last two steps, but main trick is that $\mathbf{e'M^0X = e'X=0}$)<br>
                """,
            unsafe_allow_html=True,
        )

    with st.expander("Relating two formulations of AIC (Greene pp. 47 and 561)"):
        st.markdown(
            r"""
            Not sure if this is useful, but it clarified things in my head.<br>
            
            Recall, $SSE = \mathbf{e'e}$<br>
            In the linear model with normally distributed disturbances, the maximized log likelihood is<br>
            $\text{ln} L = -\frac{n}{2} [1 + \text{ln}(2 \pi) + \text{ln}(\frac{SSE}{n})]$<br>
            Ignore the constants and notice that<br>
            $\text{ln} L \propto -\frac{n}{2} \text{ln}(\frac{SSE}{n})$<br>
            $-2 \text{ln} L \propto n \text{ln}(\frac{SSE}{n})$<br>
            $-2 \text{ln} L + 2K \propto n \text{ln}(\frac{SSE}{n}) + 2K$<br>
            $-2 \text{ln} L + 2K \propto \text{ln}(\frac{SSE}{n}) + \frac{2K}{n}$<br>
            Which we wanted to show.<br>
            Might have been enough to just state that $\text{ln} L \propto -\text{ln}(\frac{SSE}{n})$.
""",
            unsafe_allow_html=True,
        )

    # with st.expander(r"$0 \leq R^2 \leq 1$"):
    #     st.markdown("TBD")

    # with st.expander("R-sq increases with the number of regressors"):
    #     st.markdown("TBD")

pages/2_fit.py
justgri/phd_econometrics
import matplotlib.cm as cm
import plotly.express as px
import plotly.graph_objects as go
import plotly.io as pio

# Sample template from plotly documentation
# https://plotly.com/python/templates/
pio.templates["draft"] = go.layout.Template(
    layout_annotations=[
        dict(
            name="draft watermark",
            text="DRAFT",
            textangle=-30,
            opacity=0.1,
            font=dict(color="black", size=100),
            xref="paper",
            yref="paper",
            x=0.5,
            y=0.5,
            showarrow=False,
        )
    ]
)

my_font = "Sans-Serif"
my_font_color = "black"

# Personal template for streamlit charts
pio.templates["my_streamlit"] = go.layout.Template(
    layout=dict(
        # Font options - streamlit overwrites those based on config.toml file
        font={
            "family": my_font,
            "size": 12,
            "color": my_font_color,
        },
        # Title options
        title=dict(
            # text="",
            font={
                "family": my_font,
                "size": 25,
                "color": my_font_color,
            },
            x=0.00,
            xanchor="left",
            xref="paper",
            y=0.98,
            yanchor="top",
            yref="container",
        ),
        # Legend options
        showlegend=True,
        legend=dict(
            orientation="h",
            x=-0.01,
            xanchor="left",
            y=1.00,
            yanchor="bottom",
            # bgcolor="red",
            # bordercolor="#333",
            # borderwidth=2,
            entrywidth=0,
            entrywidthmode="pixels",
            font={
                "family": my_font,
                "size": 12,
                "color": my_font_color,
            },
            itemclick="toggle",
            # tracegroupgap=10,
            traceorder="normal",
            valign="middle",
        ),
        # Paper options
        margin=dict(autoexpand=True, l=70, r=40, t=25, b=40, pad=0),
        width=500,
        height=400,
        autosize=False,
        paper_bgcolor="white",  # doesn't work with streamlit
        plot_bgcolor="white",  # doesn't work with streamlit
        separators=".,",
        # Chart interaction options
        dragmode=False,  # zoom
        # Hover options
        hovermode="closest",
        hoverlabel=dict(
            font=dict(family=my_font, color=my_font_color),
            # bgcolor="rgb(255, 255, 255)",  # rgb(0, 0, 0, 0)
            bordercolor="black",
            namelength=-1,
        ),
        # Gride options
        grid=dict(columns=2, rows=2, domain=dict(x=[0, 1], y=[0, 1])),
        # Axes options
        xaxis=dict(
            automargin=False,
            autorange=True,
            categoryorder="sum descending",
            color="black",
            gridcolor="white",
            showline=True,
            linecolor="black",
            linewidth=2,
            ticks="outside",
            tickcolor="black",
            tickformat="%b %Y",
            title=dict(standoff=0, text=""),
            zerolinecolor="black",
            zerolinewidth=2,
            layer="above traces",
            showticklabels=True,
            showspikes=True,
            spikethickness=1,
        ),
        yaxis=dict(
            automargin="left+right",
            autorange=True,
            categoryorder="sum descending",
            color="black",
            gridcolor="lightgrey",
            gridwidth=0.5,
            showline=True,
            linecolor="black",
            linewidth=2,
            ticks="",
            tickformat="$,.0f",
            title=dict(
                standoff=0,
                text="",
                font=dict(family=my_font, color=my_font_color),
            ),
            separatethousands=True,
            zerolinecolor="black",
            zerolinewidth=2,
            layer="above traces",
            rangemode="tozero",
        ),
        colorway=px.colors.qualitative.G10,
    ),
    data=dict(
        scatter=[
            go.Scatter(
                line=dict(width=4),
                marker=dict(symbol="diamond", size=8),
                mode="lines",
                connectgaps=False,
                # hovertemplate="%{fullData.name}<br>%{x}:%{y}<extra></extra>",
                hovertemplate="%{x}<br>%{y}<extra></extra>",
                xhoverformat="%b %Y",
                yhoverformat="$,.0f",
            )
        ]
    ),
)

cols_g10 = px.colors.qualitative.G10
cols_set1 = px.colors.qualitative.Set1

cols_set1_plt = cm.Set1.colors

set1_red = cols_set1_plt[0]
set1_blue = cols_set1_plt[1]
set1_green = cols_set1_plt[2]
set1_purple = cols_set1_plt[3]
set1_orange = cols_set1_plt[4]
set1_yellow = cols_set1_plt[5]
set1_brown = cols_set1_plt[6]
set1_pink = cols_set1_plt[7]
set1_grey = cols_set1_plt[8]

"""
Color templates
https://plotly.com/python/discrete-color/

Personal favorites

px.colors.qualitative.G10
['#3366CC', '#DC3912', '#FF9900', '#109618', '#990099',
    '#0099C6', '#DD4477', '#66AA00', '#B82E2E', '#316395']

px.colors.qualitative.Set1
px.colors.qualitative.D3
"""

src/scripts/plot_themes.py
justgri/phd_econometrics
import streamlit as st


def local_css(file_path):
    with open(file_path) as f:
        st.markdown(f"<style>{f.read()}</style>", unsafe_allow_html=True)


def external_css(file_url):
    st.markdown(f'<link href="{file_url}" rel="stylesheet">', unsafe_allow_html=True)


def wide_col():
    return st.columns((0.02, 1, 0.02))


def narrow_col():
    return st.columns((0.35, 1, 0.35))


def narrow_col_intro():
    return st.columns((0.1, 1, 0.1))


def two_cols():
    return st.columns((0.1, 1, 0.2, 1, 0.1))

src/scripts/utils.py
justgri/phd_econometrics