Statistics and Financial Data Analysis¶

Problem Sheet 1 — Solutions¶

$\bullet$ Note this is the personal version of Wen Su rather than a university official solution.

Q1(a) — Data Preparation¶

$\bullet$ Load the data and display

In [1]:
# Required packages
import numpy as np
import pandas as pd

# Read the data
price_df = pd.read_csv("Stock_Prices.csv").set_index('date')
factor_df = pd.read_csv("Fama_French_Factors.csv").set_index('date')

# Display the data
display(price_df.head())
display(factor_df.head())
MMM ABT ADBE AMD AFL APD ALK LNT MO AEP ... WEC WFC WELL WDC WY WHR WMB XEL XRX ZION
date
1990-01-02 4.945746 1.809341 1.188340 4.1250 0.638076 5.295442 4.663832 2.243539 0.541092 5.801660 ... 2.805641 1.084990 1.308260 3.462067 3.274985 12.320812 1.465829 4.096441 10.246303 1.761991
1990-01-03 4.999504 1.815872 1.247023 4.0000 0.624499 5.295442 4.609285 2.301362 0.549003 5.779771 ... 2.750843 1.037035 1.308260 3.662765 3.202847 12.368569 1.498097 4.083797 10.335008 1.793455
1990-01-04 5.037902 1.812606 1.305707 3.9375 0.610923 5.228580 4.582011 2.289797 0.530017 5.692196 ... 2.739885 1.025046 1.330436 3.712940 3.173993 12.511831 1.452001 3.982649 10.357187 1.793455
1990-01-05 4.953426 1.793011 1.335048 3.8125 0.610923 5.295442 4.663832 2.243539 0.533181 5.560837 ... 2.696045 0.995073 1.319348 3.813290 3.159566 12.129787 1.456611 3.944721 10.379369 1.793455
1990-01-08 5.083982 1.802203 1.352692 3.8125 0.619974 5.348928 4.702105 2.197280 0.518942 5.538945 ... 2.674126 1.025046 1.319348 3.813290 3.159566 12.320812 1.475049 3.906791 10.312830 1.793455

5 rows × 236 columns

r_market rfr hml smb
date
1990-01-02 1.466 0.026 -0.07 -0.68
1990-01-03 -0.034 0.026 -0.26 0.73
1990-01-04 -0.684 0.026 -0.22 0.42
1990-01-05 -0.824 0.026 -0.21 0.77
1990-01-08 0.326 0.026 -0.26 -0.41

$\bullet$ Restrict the dataset to the period from 1 Jan 2019 to 31 Dec 2023

$\bullet$ And calculate daily return by $ r_{i,t}=\frac{P_{i,t}}{P_{i,t-1}}-1. $

$\bullet$ Excess return is $R_{i,t}=r_{i,t}-r_{f,t}$

In [2]:
price_df.index = pd.to_datetime(price_df.index)
factor_df.index = pd.to_datetime(factor_df.index)
start_date = "2019-01-01"
end_date = "2023-12-31"

return_df_5y = price_df.pct_change().loc[start_date:end_date]
factor_df_5y = factor_df.loc[start_date:end_date] / 100 # Factor Return
ereturn_df_5y = return_df_5y.sub(factor_df_5y['rfr'], axis=0) # Excess Return
factor_df_5y['mkt_excess'] = factor_df_5y['r_market'] - factor_df_5y['rfr'] # Market excess return

$\bullet$ Report summary statistics for the stock returns and factor variables. Pool all stock-day observations when summarising stock returns.

In [3]:
# Summary statistics
summary_df = pd.DataFrame({
    'Stock returns (pooled)': return_df_5y.stack().describe(),
    'Market return': factor_df_5y['r_market'].describe(),
    'Risk-free rate': factor_df_5y['rfr'].describe(),
    'SMB': factor_df_5y['smb'].describe(),
    'HML': factor_df_5y['hml'].describe()
}).T

display(summary_df.style.format('{:.2%}').format('{:.0f}', subset=['count']))
  count mean std min 25% 50% 75% max
Stock returns (pooled) 296888 0.07% 2.33% -53.86% -0.94% 0.08% 1.06% 74.59%
Market return 1258 0.07% 1.38% -11.99% -0.52% 0.10% 0.76% 9.35%
Risk-free rate 1258 0.01% 0.01% 0.00% 0.00% 0.01% 0.01% 0.02%
SMB 1258 -0.00% 0.73% -3.56% -0.45% -0.02% 0.41% 5.46%
HML 1258 -0.00% 1.13% -5.02% -0.62% -0.04% 0.61% 6.73%

$\bullet$ The resulting dataset contains daily returns for 236 stocks over 1,258 trading days, from 2 Jan 2019 to 29 Dec 2023, giving 296,888 stock-day observations. It also includes the daily market return, risk-free rate, and the SMB and HML factors. Stock return statistics are pooled across stocks and dates.

$\bullet$ The pooled stock returns have a daily mean of 0.07% and a standard deviation of 2.33%. The market return has a similar mean of 0.07% but a lower standard deviation of 1.38%, indicating greater dispersion in the pooled individual stock returns.

$\bullet$ SMB and HML have daily means close to zero, with standard deviations of 0.73% and 1.13%, respectively. HML therefore exhibits greater variability than SMB over this period. The daily risk-free rate is much smaller in both magnitude and variability than the stock returns and factor returns.

Q1(b) — CAPM and Fama–French Models for JPM¶

$\bullet$ Estimate the CAPM, where $R_{JPM,t}$ denotes JPM's excess return: $R_{JPM,t} = \alpha + \beta_M(r_{m,t}-r_{f,t}) + \varepsilon_t.$

$\bullet$ Estimate the Fama–French three-factor model: $R_{JPM,t} = \alpha + \beta_M(r_{m,t}-r_{f,t}) + \beta_S SMB_t + \beta_H HML_t + \varepsilon_t.$

In [4]:
import statsmodels.api as sm
from statsmodels.iolib.summary2 import summary_col

# HAC (Newey–West): 5 lags, Bartlett kernel
hac_kwds = {
    'maxlags': 5,
    'kernel': 'bartlett',
    'use_correction': True
}

# CAPM
X_capm = sm.add_constant(factor_df_5y[['mkt_excess']])
y = ereturn_df_5y['JPM']
capm_jpm = sm.OLS(y, X_capm).fit(
    cov_type='HAC', cov_kwds=hac_kwds, use_t=True
)

# Fama–French three-factor model
X_ff = sm.add_constant(factor_df_5y[['mkt_excess', 'smb', 'hml']])
ff_jpm = sm.OLS(y, X_ff).fit(
    cov_type='HAC', cov_kwds=hac_kwds, use_t=True
)

# Regression table
regression_df = summary_col(
    [capm_jpm, ff_jpm],
    model_names=['CAPM (HAC)', 'Fama–French (HAC)'],
    regressor_order=['const', 'mkt_excess', 'smb', 'hml'],
    stars=True,
    float_format='%.4f',
    info_dict={'N': lambda model: str(int(model.nobs))}
).tables[0]

display(regression_df.T)
const mkt_excess smb hml R-squared R-squared Adj. N
CAPM (HAC) 0.0001 (0.0004) 1.0596*** (0.0614) 0.5274 0.5270 1258
Fama–French (HAC) 0.0001 (0.0003) 1.1110*** (0.0268) -0.1463*** (0.0560) 0.8594*** (0.0490) 0.7617 0.7611 1258

$\bullet$ A one-percentage-point increase in the market excess return is associated with an increase of about 1.06 percentage points in JPM's excess return under CAPM, or 1.11 under Fama--French holding SMB and HML fixed. Both market coefficients are significant at the 1\% level.

$\bullet$ The SMB loading is negative ($-0.1463$) and the HML loading is positive ($0.8594$). Both are significant at the 1\% level, controlling for the other factors.

$\bullet$ Neither model has a significant alpha at the 5\% level using HAC standard errors ($p=0.8396$ for CAPM and $p=0.8109$ for Fama--French).

$\bullet$ The additional factors improve the in-sample fit: $R^2$ increases from $0.5274$ to $0.7617$, and adjusted $R^2$ from $0.5270$ to $0.7611$.

$\bullet$ A zero alpha means that JPM's expected excess return is explained by its factor exposures. However, failure to reject $\alpha=0$ does not establish that the true alpha is zero.

Q1(c) — Market Betas and Average Excess Returns¶

$\bullet$ Estimate the CAPM for each stock and calculate its average daily excess return. Also save the alpha estimates and their $p$-values for part (d).

In [5]:
# CAPM for all stocks
capm_results = []
for stock in ereturn_df_5y.columns:
    model = sm.OLS(ereturn_df_5y[stock], X_capm).fit(
        cov_type='HAC', cov_kwds=hac_kwds, use_t=True
    )
    capm_results.append([
        model.params['mkt_excess'],
        model.params['const'],
        model.pvalues['const']
    ])

capm_df = pd.DataFrame(
    capm_results,
    index=ereturn_df_5y.columns,
    columns=['beta', 'alpha', 'alpha_pvalue']
)
capm_df['mean_ereturn'] = ereturn_df_5y.mean()

display(capm_df.head().round(6))
beta alpha alpha_pvalue mean_ereturn
MMM 0.776525 -0.000620 0.086854 -0.000166
ABT 0.780011 0.000014 0.968532 0.000469
ADBE 1.246271 0.000244 0.586627 0.000972
AMD 1.591567 0.001220 0.081185 0.002150
AFL 0.943728 0.000144 0.729462 0.000695

$\bullet$ Plot average excess returns against estimated market betas. Add a fitted line to show the cross-sectional relationship.

In [6]:
import matplotlib.pyplot as plt
from matplotlib.ticker import PercentFormatter

# Annualize average daily excess returns
beta = capm_df['beta']
mean_ereturn = capm_df['mean_ereturn'] * 252

# Fit a straight line
slope, intercept = np.polyfit(beta, mean_ereturn, 1)
beta_grid = np.linspace(beta.min(), beta.max(), 100)
correlation = beta.corr(mean_ereturn)

# Plot settings
with plt.rc_context({
    'font.family': 'serif',
    'font.serif': ['TeX Gyre Pagella', 'DejaVu Serif'],
    'font.size': 11,
    'axes.labelsize': 12,
    'axes.titlesize': 13,
    'pdf.fonttype': 42
}):
    fig, ax = plt.subplots(figsize=(7.5, 4.8), dpi=180)

    # Stocks and fitted line
    ax.scatter(beta, mean_ereturn, s=28, color='0.55',
               alpha=0.60, linewidths=0, label='Stocks', zorder=3)
    ax.plot(beta_grid, intercept + slope * beta_grid,
            color='tab:blue', linewidth=2, label='Fitted line', zorder=4)

    # Zero-return reference and light grid
    ax.axhline(0, color='0.45', linestyle='--', linewidth=0.8, zorder=2)
    ax.set_axisbelow(True)
    ax.grid(axis='y', color='0.88', linewidth=0.6, alpha=0.7)

    # Labels and axes
    ax.set_xlabel(r'Estimated market beta, $\hat{\beta}_i$')
    ax.set_ylabel('Annualized average excess return')
    ax.set_title('Market beta and average excess return', pad=12)
    ax.yaxis.set_major_formatter(PercentFormatter(xmax=1, decimals=0))
    ax.tick_params(direction='out', length=4, width=0.8)
    ax.margins(x=0.05, y=0.10)

    for spine in ['top', 'right']:
        ax.spines[spine].set_visible(False)

    # Legend and correlation
    ax.legend(loc='upper left', frameon=False, fontsize=10)
    ax.text(0.97, 0.05, f'Correlation = {correlation:.3f}',
            transform=ax.transAxes, ha='right', va='bottom',
            fontsize=10, color='0.35')

    plt.tight_layout()
    plt.show()

$\bullet$ The relationship is positive: stocks with higher estimated betas tend to have higher average excess returns. The correlation is $0.375$, but there is substantial dispersion around the fitted line.

Q1(d) — Comparing Alphas¶

$\bullet$ Estimate the Fama–French model for each stock and collect the alpha estimates and their $p$-values.

In [7]:
# Fama–French model for all stocks
ff_results = []
for stock in ereturn_df_5y.columns:
    model = sm.OLS(ereturn_df_5y[stock], X_ff).fit(
        cov_type='HAC', cov_kwds=hac_kwds, use_t=True
    )
    ff_results.append([
        model.params['const'],
        model.pvalues['const']
    ])

ff_df = pd.DataFrame(
    ff_results,
    index=ereturn_df_5y.columns,
    columns=['alpha', 'alpha_pvalue']
)

# Alpha summary statistics
alpha_df = pd.DataFrame({
    'CAPM': capm_df['alpha'],
    'Fama–French': ff_df['alpha']
})
alpha_summary = alpha_df.agg(['mean', 'median', 'std', 'min', 'max'])
alpha_summary.loc['Mean absolute alpha'] = alpha_df.abs().mean()

# Unadjusted individual two-sided HAC tests at the 5% level
significant_capm = capm_df['alpha_pvalue'] < 0.05
significant_ff = ff_df['alpha_pvalue'] < 0.05

alpha_summary.loc['Significant alphas'] = [
    significant_capm.sum(), significant_ff.sum()
]
alpha_summary.loc['Proportion significant'] = [
    significant_capm.mean(), significant_ff.mean()
]

display(alpha_summary.round(6).T)
mean median std min max Mean absolute alpha Significant alphas Proportion significant
CAPM 0.000041 0.000041 0.000374 -0.001260 0.001220 0.000282 5.0 0.021186
Fama–French 0.000048 0.000046 0.000367 -0.001213 0.001274 0.000277 6.0 0.025424

$\bullet$ Using HAC standard errors, CAPM has 5 significant alphas out of 236 stocks ($2.12\%$), while Fama--French has 6 ($2.54\%$). Thus, the extra factors do not reduce the proportion of significant alphas in this sample. Positive and negative alphas can cancel in the mean. The mean absolute alpha decreases only slightly, from $0.000282$ under CAPM to $0.000277$ under Fama--French.

$\bullet$ These are unadjusted individual tests at the $5\%$ level. Testing many stocks creates a multiple-testing issue, so some significant alphas may be false positives. A zero alpha means that the asset's expected excess return is explained by its factor exposures; failure to reject $\alpha=0$ does not establish that the true alpha is zero.

Q2(a) — Ordinary Least Squares¶

$\bullet$ Load the data into df_WLS and regress.

In [8]:
# Required packages
import numpy as np
import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt

# Read the data
df_WLS = pd.read_csv("WLS_Data.txt", index_col=0)

# OLS regression
x = df_WLS['x']
y = df_WLS['y_data']
X = sm.add_constant(df_WLS[['x']])
ols_model = sm.OLS(y, X).fit()
In [9]:
# Plot settings
plot_style = {
    'font.family': 'serif',
    'font.serif': ['TeX Gyre Pagella', 'DejaVu Serif'],
    'font.size': 11,
    'axes.labelsize': 12,
    'axes.titlesize': 13,
    'pdf.fonttype': 42
}

with plt.rc_context(plot_style):
    fig, axes = plt.subplots(1, 2, figsize=(12, 4.5), dpi=180)

    # (a) Original data and OLS fitted values
    ax = axes[0]
    ax.scatter(x, y, s=28, color='0.55', alpha=0.75,
               linewidths=0, label='Data', zorder=3)
    ax.plot(x, ols_model.fittedvalues, color='tab:blue',
            linewidth=2, label='OLS fitted values', zorder=4)

    # Dashed lines showing residuals
    ax.vlines(x, y, ols_model.fittedvalues,
              colors='0.65', linestyles='--',
              linewidths=0.8, alpha=0.5, zorder=2)

    ax.set_xlabel('$x$')
    ax.set_ylabel('$y$')
    ax.set_title('(a) OLS fitted values and original data', pad=12)
    ax.legend(frameon=False)

    # (b) Residuals against x
    ax = axes[1]
    ax.scatter(x, ols_model.resid, s=28, color='0.55',
               alpha=0.75, linewidths=0, zorder=3)
    ax.axhline(0, color='tab:blue', linewidth=1.8, zorder=2)

    ax.set_xlabel('$x$')
    ax.set_ylabel('OLS residuals')
    ax.set_title('(b) Residuals against $x$', pad=12)

    # Shared formatting
    for ax in axes:
        ax.set_axisbelow(True)
        ax.grid(axis='y', color='0.88', linewidth=0.6, alpha=0.7)
        ax.tick_params(direction='out', length=4, width=0.8)
        for spine in ['top', 'right']:
            ax.spines[spine].set_visible(False)

    plt.tight_layout(w_pad=2.5)
    plt.show()

$\bullet$ The residuals are much more dispersed for the larger values of x, particularly among the last 20 observations.

$\bullet$ This suggests heteroskedasticity rather than constant error variance, i.e., $\mathrm{Var}(\varepsilon_i|\boldsymbol X) \neq$ constant.

Q2(b) — Weighted Least Squares¶

$\bullet$ Let $w_i$ denote df_WLS['weights'].

$\bullet$ The question states that $w_i^2 \propto \operatorname{Var}(\varepsilon_i| X)$, we assume $\operatorname{Var}(\varepsilon_i| X)=\sigma^2w_i^2$.

$\bullet$ Therefore, WLS minimises

$ \sum_{i=1}^{n}(\operatorname{Var}(\varepsilon_i| X))^{-1}{(y_i-\beta_0-\beta_1x_i)^2} \propto \sum_{i=1}^{n}w_i^{-2}{(y_i-\beta_0-\beta_1x_i)^2} $

$\bullet$ In statsmodels, use weights=1 / w**2 to implement this transformation.

$\bullet$ Compare the coefficient estimates and their estimated standard deviations (standard errors) with OLS.

In [10]:
# WLS: inverse-variance weights
w = df_WLS['weights']
wls_model = sm.WLS(y, X, weights=1 / w**2).fit()

# Coefficient estimates and standard errors
coefficient_df = pd.DataFrame({
    'OLS estimate': ols_model.params,
    'OLS std. error': ols_model.bse,
    'WLS estimate': wls_model.params,
    'WLS std. error': wls_model.bse
})

display(coefficient_df.style.format('{:.4f}'))
  OLS estimate OLS std. error WLS estimate WLS std. error
const 5.2426 0.2707 5.2469 0.1426
x 0.4349 0.0233 0.4466 0.0180

$\bullet$ The estimated intercept is $5.2426$ under OLS and $5.2469$ under WLS; the slope changes from $0.4349$ to $0.4466$. The fitted relationships are therefore similar.

$\bullet$ WLS reports smaller standard errors for both coefficients. The OLS standard errors shown here use the usual constant-variance formula, which is not appropriate under heteroskedasticity.

Q2(c) — Prediction Intervals¶

$\bullet$ The pointwise 95% prediction interval for an observation at $\mathbf{x}_i=(1,x_i)^\top$, with relative error standard deviation $w_i$, is $$ \widehat{y}_i \pm t_{n-2,\,0.975}\,\widehat{\sigma} \sqrt{w_i^2+\mathbf{x}_i^\top(X^\top W X)^{-1}\mathbf{x}_i}, $$ where $W=\operatorname{diag}(w_1^{-2},\ldots,w_n^{-2})$ and $\widehat{\sigma}^2 =\frac{1}{n-2}\sum_{i=1}^{n}\frac{\widehat{\varepsilon}_i^2}{w_i^2}$.

In [11]:
split_x = x[w != 1].min()
x_grid = np.r_[np.linspace(x.min(), split_x, 250),
               np.linspace(split_x, x.max(), 250)]
w_grid = np.r_[np.ones(250), np.full(250, w[w != 1].iloc[0])]
X_grid = sm.add_constant(pd.DataFrame({'x': x_grid}))

ols_prediction = ols_model.get_prediction(
    X_grid
).summary_frame(alpha=0.05)

wls_prediction = wls_model.get_prediction(
    X_grid, weights=1 / w_grid**2
).summary_frame(alpha=0.05)

with plt.rc_context(plot_style):
    fig, ax = plt.subplots(figsize=(7.5, 4.8), dpi=180)

    ax.scatter(x, y, s=28, color='0.55', alpha=0.75,
               linewidths=0, label='Data', zorder=3)

    for name, prediction, color in [
        ('OLS', ols_prediction, 'tab:orange'),
        ('WLS', wls_prediction, 'tab:blue')
    ]:
        ax.plot(x_grid, prediction['mean'], color=color,
                linewidth=2, label=f'{name} fitted values')
        ax.fill_between(
            x_grid,
            prediction['obs_ci_lower'].to_numpy(),
            prediction['obs_ci_upper'].to_numpy(),
            color=color, alpha=0.15,
            label=f'{name} 95% prediction interval'
        )

    ax.set_xlabel('$x$')
    ax.set_ylabel('$y$')
    ax.set_title('OLS and WLS prediction intervals')
    ax.set_axisbelow(True)
    ax.grid(axis='y', color='0.88', linewidth=0.6)
    for spine in ['top', 'right']:
        ax.spines[spine].set_visible(False)
    ax.legend(frameon=False, fontsize=9)
    plt.tight_layout()
    plt.show()

$\bullet$ OLS and WLS give similar fitted lines. The OLS prediction intervals have roughly constant width, whereas the WLS intervals are narrower for the first 30 observations and wider for the last 20, reflecting the larger error variance in the second group.

$\bullet$ The usual OLS intervals assume constant error variance, so their nominal 95% coverage is not justified here. WLS accounts for the given variance differences.

Q3(a) — Mean Squared Error¶

$\bullet$ Consider sample variance $S^2=\frac{1}{n-1}\sum_{i=1}^n(X_i-\bar X)^2$, we have $ \theta=\alpha\sum_{i=1}^n(X_i-\bar X)^2 =\alpha(n-1)S^2.$

$\bullet$ Since $\mathbb{E}(S^2) = \sigma^2$, one has $ \mathbb{E}[\theta] =\alpha(n-1)\sigma^2. $

$\bullet$ The bias is therefore $ \operatorname{Bias}(\theta) =\mathbb{E}[\theta]-\sigma^2 =[\alpha(n-1)-1]\sigma^2. $

$\bullet$ Using the variance formula given in the question, $ \operatorname{Var}(\theta) =\alpha^2(n-1)^2\operatorname{Var}(S^2) =\alpha^2(n-1)^2 \left(\frac{\gamma}{n}+\frac{2}{n-1}\right)\sigma^4. $

$\bullet$ Combining the variance and squared bias gives $$ \boxed{ \operatorname{MSE}(\theta) =\sigma^4\left[ \alpha^2(n-1)^2\left(\frac{\gamma}{n}+\frac{2}{n-1}\right) +[\alpha(n-1)-1]^2 \right]. } $$

Q3(b) — The MSE-Minimising Choice of $\alpha$¶

$\bullet$ Differentiate the MSE with respect to $\alpha$ and set the derivative equal to zero: $$ \frac{\partial\operatorname{MSE}(\theta)}{\partial\alpha} =2\sigma^4(n-1)\left[ \alpha(n-1)\left(1+\frac{\gamma}{n}+\frac{2}{n-1}\right)-1 \right]=0. $$ Hence, $ \boxed{ \alpha^* =\frac{1}{(n-1)(1+\gamma/n)+2} =\frac{1}{n+1+\frac{n-1}{n}\gamma}. } $

$\bullet$ This is the unique minimum because $ \frac{\partial^2\operatorname{MSE}(\theta)}{\partial\alpha^2} =2\sigma^4(n-1)^2\left(1+\frac{\gamma}{n}+\frac{2}{n-1}\right)>0. $ The optimal $\alpha$ depends on $n$ and $\gamma$, but not on $\sigma^2$.

Q3(c) — The Normal Case¶

$\bullet$ In the normal case, $\gamma=0$, so $\alpha^*=1/(n+1)$. The minimum-MSE estimator is $$ \theta^*=\frac{1}{n+1}\sum_{i=1}^n(X_i-\bar X)^2 =\frac{n-1}{n+1}S^2. $$ Its bias is $ \boxed{ \operatorname{Bias}(\theta^*) =\left(\frac{n-1}{n+1}-1\right)\sigma^2 =-\frac{2\sigma^2}{n+1}. } $

$\bullet$ The variance and squared bias give $ \operatorname{MSE}(\theta^*) =\frac{2(n-1)\sigma^4}{(n+1)^2} +\frac{4\sigma^4}{(n+1)^2} =\frac{2\sigma^4}{n+1}. $

$\bullet$ For the unbiased estimator $S^2$, corresponding to $\alpha=1/(n-1)$, $ \operatorname{MSE}(S^2)=\operatorname{Var}(S^2) =\frac{2\sigma^4}{n-1}. $

$\bullet$ The reduction in MSE relative to the unbiased estimator is $$ \boxed{ \operatorname{MSE}(S^2)-\operatorname{MSE}(\theta^*) =\frac{2\sigma^4}{n-1}-\frac{2\sigma^4}{n+1} =\frac{4\sigma^4}{n^2-1}>0. } $$ Thus, use of MSE tends to be forgiving of underestimation.

Q3(d) — The Bias–Variance Trade-off¶

$\bullet$ For $n=20$, $\sigma^2=1$ and $\gamma=0$, $$ \operatorname{Bias}(\theta)^2=(19\alpha-1)^2, \qquad \operatorname{Var}(\theta)=38\alpha^2, \qquad \operatorname{MSE}(\theta)=(19\alpha-1)^2+38\alpha^2. $$

$\bullet$ Plot the three quantities against $\alpha$, marking the MSE-minimising value $\alpha^*=1/21$ and the unbiased choice $\alpha_{\mathrm{U}}=1/19$.

In [12]:
# Required packages
import numpy as np
import matplotlib.pyplot as plt

# Squared bias, variance and MSE
n = 20
sigma2 = 1
gamma = 0
alpha_grid = np.linspace(0.02, 0.08, 500)
bias_squared = (alpha_grid * (n - 1) - 1)**2 * sigma2**2
variance = alpha_grid**2 * (n - 1)**2 * (gamma / n + 2 / (n - 1)) * sigma2**2
mse = bias_squared + variance

# Optimal and unbiased choices
alpha_opt = 1 / (n + 1 + (n - 1) * gamma / n)
alpha_unbiased = 1 / (n - 1)
mse_opt = 2 * sigma2**2 / (n + 1)
mse_unbiased = 2 * sigma2**2 / (n - 1)

# Plot settings
plot_style = {
    'font.family': 'serif',
    'font.serif': ['TeX Gyre Pagella', 'DejaVu Serif'],
    'font.size': 11,
    'axes.labelsize': 12,
    'axes.titlesize': 13,
    'pdf.fonttype': 42
}

with plt.rc_context(plot_style):
    fig, ax = plt.subplots(figsize=(7.5, 4.8), dpi=180)
    ax.plot(alpha_grid, bias_squared, color='0.45', linestyle='--',
            linewidth=1.8, label='Squared bias')
    ax.plot(alpha_grid, variance, color='0.65', linestyle=':',
            linewidth=2, label='Variance')
    ax.plot(alpha_grid, mse, color='tab:blue', linewidth=2,
            label='MSE')

    # Mark the optimal and unbiased choices
    ax.axvline(alpha_opt, ymax=0.82, color='tab:blue', linestyle='--',
               linewidth=0.8, alpha=0.55)
    ax.axvline(alpha_unbiased, ymax=0.82, color='0.45', linestyle=':',
               linewidth=0.9, alpha=0.8)
    ax.scatter(alpha_opt, mse_opt, color='tab:blue', s=30, zorder=4)
    ax.scatter(alpha_unbiased, mse_unbiased, color='0.45',
               s=30, marker='D', zorder=4)
    ax.text(alpha_opt - 0.0006, 0.85, r'$\alpha^*=1/21$',
            transform=ax.get_xaxis_transform(), ha='right',
            fontsize=10, color='tab:blue')
    ax.text(alpha_unbiased + 0.0006, 0.85, r'$\alpha_{\mathrm{U}}=1/19$',
            transform=ax.get_xaxis_transform(), ha='left',
            fontsize=10, color='0.35')

    ax.set_xlabel(r'$\alpha$')
    ax.set_ylabel('Squared bias, variance and MSE')
    ax.set_title('Bias–variance trade-off', pad=12)
    ax.set_ylim(bottom=0)
    ax.set_axisbelow(True)
    ax.grid(axis='y', color='0.88', linewidth=0.6, alpha=0.7)
    ax.tick_params(direction='out', length=4, width=0.8)
    for spine in ['top', 'right']:
        ax.spines[spine].set_visible(False)
    ax.legend(loc='upper center', ncol=3, frameon=False, fontsize=10)
    plt.tight_layout()
    plt.show()

$\bullet$ Squared bias is zero at $\alpha_{\mathrm{U}}=1/19$. Reducing $\alpha$ to $\alpha^*=1/21$ increases squared bias from $0$ to $4/441\approx0.0091$, but reduces the variance from $2/19\approx0.1053$ to $38/441\approx0.0862$. The reduction in variance is larger, so the MSE falls from $0.1053$ to $0.0952$.

$\bullet$ Reducing $\alpha$ further lowers the variance but increases squared bias enough to raise the MSE. The MSE is therefore minimised at a biased estimator rather than at the unbiased choice.