Lab 4 Homework: Volatility Modelling with GARCH

Understanding and Estimating Time-Varying Volatility

NoteExpected Time
  • Core homework: 90–120 minutes
  • Optional extensions: +30 minutes

Complete this homework before the in-class Bloomberg session.

Open in Colab
NoteIn-Class Bloomberg Session

After completing this homework, bring your results to the Lab 4 Bloomberg in-class session, where you will extract VIX and SPY data directly from the Bloomberg Terminal and compare implied vs realised volatility.

1 Learning Objectives

By the end of this homework, you should be able to:

  1. Identify volatility clustering in financial return data
  2. Estimate GARCH(1,1) models using the arch package
  3. Interpret GARCH parameters and their implications
  4. Compare symmetric vs asymmetric volatility models
  5. Evaluate model fit using diagnostic tools

2 Setup

Show code 
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats

# Install arch if needed (Colab)
try:
    from arch import arch_model
except ImportError:
    import subprocess
    subprocess.check_call(['pip', 'install', 'arch'])
    from arch import arch_model

np.random.seed(42)

3 Part 1: Understanding Volatility Clustering

3.1 Task 1.1: Simulate Returns with GARCH Dynamics

Before working with real data, let’s understand GARCH by simulating it.

Show code 
def simulate_garch(n, omega, alpha, beta, df=5):
    """
    Simulate returns from a GARCH(1,1) process with t-distributed innovations.
    
    Parameters
    ----------
    n : int
        Number of observations
    omega : float
        Constant term in variance equation
    alpha : float
        ARCH parameter (reaction to shocks)
    beta : float
        GARCH parameter (persistence)
    df : float
        Degrees of freedom for t-distribution
        
    Returns
    -------
    returns : array
        Simulated returns
    sigma2 : array
        True conditional variance
    """
    sigma2 = np.zeros(n)
    returns = np.zeros(n)
    
    # Unconditional variance as starting value
    sigma2[0] = omega / (1 - alpha - beta)
    
    for t in range(1, n):
        sigma2[t] = omega + alpha * returns[t-1]**2 + beta * sigma2[t-1]
        returns[t] = np.sqrt(sigma2[t]) * np.random.standard_t(df)
    
    return returns, sigma2

# Simulate 1000 daily returns
n = 1000
omega, alpha, beta = 0.00001, 0.08, 0.90
returns_sim, sigma2_true = simulate_garch(n, omega, alpha, beta)

3.2 Task 1.2: Visualise Volatility Clustering

Your task: Create a figure with three panels showing: 1. The simulated returns 2. Absolute returns (proxy for volatility) 3. True conditional standard deviation

Show code 
fig, axes = plt.subplots(3, 1, figsize=(12, 8), sharex=True)

# Panel 1: Returns
axes[0].plot(returns_sim, linewidth=0.5, color='steelblue')
axes[0].axhline(0, color='gray', linestyle='--', linewidth=0.5)
axes[0].set_ylabel('Return')
axes[0].set_title('Simulated GARCH(1,1) Returns')

# Panel 2: Absolute returns
axes[1].plot(np.abs(returns_sim), linewidth=0.5, color='coral')
axes[1].set_ylabel('|Return|')
axes[1].set_title('Absolute Returns (Volatility Proxy)')

# Panel 3: True conditional volatility
axes[2].plot(np.sqrt(sigma2_true), linewidth=1, color='darkgreen')
axes[2].set_ylabel('σ')
axes[2].set_xlabel('Time')
axes[2].set_title('True Conditional Standard Deviation')

plt.tight_layout()
plt.show()
TipReflection 1.1

Look at the three panels. Can you identify periods of high and low volatility? How do the absolute returns compare to the true conditional volatility?

3.3 Task 1.3: ACF Analysis

Your task: Compare the ACF of returns vs squared returns.

Show code 
from statsmodels.graphics.tsaplots import plot_acf

fig, axes = plt.subplots(1, 2, figsize=(12, 4))

# ACF of returns
plot_acf(returns_sim, ax=axes[0], lags=30, title='ACF: Returns')

# ACF of squared returns  
plot_acf(returns_sim**2, ax=axes[1], lags=30, title='ACF: Squared Returns')

plt.tight_layout()
plt.show()
TipReflection 1.2

What pattern do you see in the two ACF plots? Why does this pattern support the use of GARCH models?

4 Part 2: Estimating GARCH Models

4.1 Task 2.1: Fit GARCH(1,1)

Now let’s estimate the parameters from the simulated data and see how well we recover the true values.

Show code 
# Scale returns to percentages for numerical stability
returns_pct = returns_sim * 100

# Fit GARCH(1,1)
model = arch_model(returns_pct, vol='Garch', p=1, q=1, mean='Constant')
result = model.fit(disp='off')

print("=== GARCH(1,1) Estimation Results ===")
print(result.summary().tables[1])

4.2 Task 2.2: Compare Estimated vs True Parameters

Show code 
# Extract estimated parameters (note: scaled by 100^2 for variance)
omega_est = result.params['omega'] / 10000
alpha_est = result.params['alpha[1]']
beta_est = result.params['beta[1]']

print("\n=== Parameter Comparison ===")
print(f"{'Parameter':<12} {'True':>10} {'Estimated':>12} {'Error (%)':>12}")
print("-" * 48)
print(f"{'omega':<12} {omega:.6f} {omega_est:>12.6f} {(omega_est-omega)/omega*100:>11.1f}%")
print(f"{'alpha':<12} {alpha:>10.2f} {alpha_est:>12.4f} {(alpha_est-alpha)/alpha*100:>11.1f}%")
print(f"{'beta':<12} {beta:>10.2f} {beta_est:>12.4f} {(beta_est-beta)/beta*100:>11.1f}%")
print(f"{'alpha+beta':<12} {alpha+beta:>10.2f} {alpha_est+beta_est:>12.4f}")
TipReflection 2.1

How well did the estimation recover the true parameters? What does the persistence (α + β) tell us about this process?

4.3 Task 2.3: Compare Estimated vs True Volatility

Show code 
# Get conditional volatility from the model
cond_vol_est = result.conditional_volatility / 100  # Scale back

fig, ax = plt.subplots(figsize=(12, 5))
ax.plot(np.sqrt(sigma2_true), label='True σ', color='darkgreen', alpha=0.7)
ax.plot(cond_vol_est, label='Estimated σ', color='coral', linewidth=1.5)
ax.legend()
ax.set_xlabel('Time')
ax.set_ylabel('Conditional Standard Deviation')
ax.set_title('True vs Estimated Conditional Volatility')
plt.tight_layout()
plt.show()

# Calculate correlation
corr = np.corrcoef(np.sqrt(sigma2_true[1:]), cond_vol_est[1:])[0,1]
print(f"Correlation between true and estimated volatility: {corr:.4f}")

5 Part 3: Asymmetric Models (GJR-GARCH)

5.1 Task 3.1: Simulate Asymmetric Volatility

Real markets exhibit the leverage effect: negative returns increase volatility more than positive returns.

Show code 
def simulate_gjr_garch(n, omega, alpha, beta, gamma, df=5):
    """
    Simulate returns from a GJR-GARCH(1,1) process.
    
    gamma > 0 implies leverage effect (negative shocks have larger impact)
    """
    sigma2 = np.zeros(n)
    returns = np.zeros(n)
    
    sigma2[0] = omega / (1 - alpha - beta - gamma/2)
    
    for t in range(1, n):
        indicator = 1 if returns[t-1] < 0 else 0
        sigma2[t] = omega + (alpha + gamma * indicator) * returns[t-1]**2 + beta * sigma2[t-1]
        returns[t] = np.sqrt(sigma2[t]) * np.random.standard_t(df)
    
    return returns, sigma2

# Simulate with leverage effect
np.random.seed(123)
gamma = 0.10  # Asymmetry parameter
returns_gjr, sigma2_gjr = simulate_gjr_garch(n, omega, alpha, beta, gamma)

5.2 Task 3.2: Fit Both Models and Compare

Show code 
returns_gjr_pct = returns_gjr * 100

# Fit standard GARCH
model_garch = arch_model(returns_gjr_pct, vol='Garch', p=1, q=1, mean='Constant')
result_garch = model_garch.fit(disp='off')

# Fit GJR-GARCH
model_gjr = arch_model(returns_gjr_pct, vol='Garch', p=1, o=1, q=1, mean='Constant')
result_gjr = model_gjr.fit(disp='off')

print("=== Model Comparison ===")
print(f"{'Model':<15} {'Log-Likelihood':>15} {'AIC':>12} {'BIC':>12}")
print("-" * 56)
print(f"{'GARCH(1,1)':<15} {result_garch.loglikelihood:>15.2f} {result_garch.aic:>12.2f} {result_garch.bic:>12.2f}")
print(f"{'GJR-GARCH(1,1)':<15} {result_gjr.loglikelihood:>15.2f} {result_gjr.aic:>12.2f} {result_gjr.bic:>12.2f}")

5.3 Task 3.3: Examine the Asymmetry Parameter

Show code 
print("\n=== GJR-GARCH Parameter Estimates ===")
print(result_gjr.summary().tables[1])

print(f"\nTrue gamma: {gamma:.2f}")
print(f"Estimated gamma: {result_gjr.params['gamma[1]']:.4f}")
TipReflection 3.1
  • Does the GJR model find evidence of the leverage effect (γ > 0)?
  • Which model has the better AIC/BIC? What does this tell us?
  • Why might the leverage effect matter for risk management?

6 Part 4: Working with Real Data (Bloomberg)

Parts 1–3 used simulated data so you could see exactly how GARCH works when you know the truth. Now apply the same methods to real SPY data from our shared Bloomberg database.

6.1 Task 4.1: Load Bloomberg Data

Show code 
# Load SPY from shared Bloomberg database
# Falls back automatically if no local file is available
url = "https://raw.githubusercontent.com/quinfer/fin510-colab-notebooks/main/resources/bloomberg_database.parquet"
full = pd.read_parquet(url)

spy_prices = (
    full[full['ticker'] == 'SPY'][['date', 'PX_LAST']]
    .rename(columns={'date': 'Date', 'PX_LAST': 'SPY'})
    .set_index('Date')
    .sort_index()
)

# Calculate percentage returns
returns_spy = spy_prices['SPY'].pct_change().dropna() * 100

print(f"Loaded {len(returns_spy)} daily SPY returns")
print(f"\nSummary Statistics:")
print(f"  Mean:     {returns_spy.mean():.3f}%")
print(f"  Std Dev:  {returns_spy.std():.3f}%")
print(f"  Skewness: {stats.skew(returns_spy):.3f}")
print(f"  Kurtosis: {stats.kurtosis(returns_spy):.3f}")

6.2 Task 4.2: Fit Models to Real Data

Show code 
# Fit GARCH(1,1)
model_spy_garch = arch_model(returns_spy, vol='Garch', p=1, q=1, mean='Constant')
result_spy_garch = model_spy_garch.fit(disp='off')

# Fit GJR-GARCH(1,1)
model_spy_gjr = arch_model(returns_spy, vol='Garch', p=1, o=1, q=1, mean='Constant')
result_spy_gjr = model_spy_gjr.fit(disp='off')

print("=== SPY Volatility Model Comparison ===")
print(result_spy_gjr.summary())

6.3 Task 4.3: Visualise Conditional Volatility

Show code 
fig, axes = plt.subplots(2, 1, figsize=(12, 8), sharex=True)

axes[0].plot(returns_spy.index, returns_spy.values, linewidth=0.5, color='steelblue')
axes[0].axhline(0, color='gray', linestyle='--', linewidth=0.5)
axes[0].set_ylabel('Return (%)')
axes[0].set_title('SPY Daily Returns (Bloomberg)')

vol_garch = result_spy_garch.conditional_volatility
vol_gjr = result_spy_gjr.conditional_volatility

axes[1].plot(vol_garch.index, vol_garch.values, label='GARCH(1,1)', alpha=0.7)
axes[1].plot(vol_gjr.index, vol_gjr.values, label='GJR-GARCH(1,1)', alpha=0.7)
axes[1].set_ylabel('Conditional Volatility (%)')
axes[1].set_xlabel('Date')
axes[1].set_title('Estimated Conditional Volatility')
axes[1].legend()

plt.tight_layout()
plt.show()
TipReflection 4.1

Compare your results on real data to the simulated results from Parts 1–3:

  • Are the estimated parameters (α, β) similar to the simulated ones? What does that tell you?
  • Is the leverage effect (γ) present in real SPY data?
  • Can you identify the COVID-19 shock (March 2020) in the conditional volatility plot?

7 Part 5: Diagnostic Checks

7.1 Task 5.1: Standardised Residuals

Good volatility models should produce standardised residuals that are approximately i.i.d.

Show code 
# Using simulated GJR data for demonstration
std_resid = result_gjr.std_resid

fig, axes = plt.subplots(2, 2, figsize=(12, 8))

# Histogram
axes[0, 0].hist(std_resid, bins=50, density=True, alpha=0.7, color='steelblue')
x = np.linspace(-4, 4, 100)
axes[0, 0].plot(x, stats.norm.pdf(x), 'r-', lw=2, label='N(0,1)')
axes[0, 0].plot(x, stats.t.pdf(x, df=5), 'g--', lw=2, label='t(5)')
axes[0, 0].legend()
axes[0, 0].set_title('Distribution of Standardised Residuals')

# Q-Q plot
stats.probplot(std_resid, dist="norm", plot=axes[0, 1])
axes[0, 1].set_title('Q-Q Plot (vs Normal)')

# ACF of residuals
plot_acf(std_resid, ax=axes[1, 0], lags=20, title='ACF: Standardised Residuals')

# ACF of squared residuals
plot_acf(std_resid**2, ax=axes[1, 1], lags=20, title='ACF: Squared Standardised Residuals')

plt.tight_layout()
plt.show()
TipReflection 5.1

Examine the diagnostic plots: - Do the standardised residuals look i.i.d.? - Is there remaining autocorrelation in squared residuals? - Does the distribution look normal or heavy-tailed?

8 Homework Deliverables

Complete the following before the in-class Bloomberg session:

  1. Summary Table: Fill in the table below with your results
Metric Simulated GARCH Simulated GJR Real Data (SPY)
α (ARCH)
β (GARCH)
γ (if GJR)
Persistence (α+β)
AIC
  1. Written Reflection (150-200 words): Explain in your own words:
    • Why do we model volatility as time-varying rather than constant?
    • What does high persistence (α + β close to 1) imply for forecasting?
    • When would you choose GJR over standard GARCH?
  2. Prepare for In-Class: The Bloomberg session will compare VIX (implied volatility) to realised volatility. Think about:
    • What is the VIX measuring?
    • Why might VIX differ from GARCH-estimated volatility?

9 Optional extensions

Advanced students should also complete:

  1. Model Selection: Fit GARCH(1,1), GARCH(2,1), and EGARCH(1,1) to the SPY data. Use AIC/BIC to select the best model. Is the ranking consistent?

  2. Forecasting: Use the fitted model to produce 10-day ahead volatility forecasts. Compare these to realised volatility over the next 10 days.

  3. News Impact Curves: Plot the news impact curve for both GARCH and GJR-GARCH. Quantify the asymmetry.

10 References

  • Tsay (2010) Chapter 3 on Conditional Heteroscedastic Models
  • Brooks (2019) Chapter 9 on Modelling and Forecasting Volatility

References

Brooks, Chris. 2019. Introductory Econometrics for Finance. 4th ed. Cambridge, UK: Cambridge University Press.
Tsay, Ruey S. 2010. Analysis of Financial Time Series. 3rd ed. Hoboken, NJ: Wiley.