Introduction to Statistics¶

Problem Sheet — Solutions¶

Oxford MSc Mathematical and Computational Finance

Part A¶

Q1 — Download and Prepare Adjusted Stock Prices¶

$\bullet$ Download daily adjusted closing prices for Tesla (TSLA) and Goldman Sachs (GS) from 1 Jan 2013 to 14 Aug 2020, inclusive, and combine them into a single DataFrame.

$\bullet$ Notes. (a) The problem sheet incorrectly lists Tesla's ticker as TESL; the correct ticker is TSLA. (b) We use yfinance instead of pandas_datareader, as the latter's Yahoo Finance interface is no longer reliable.

$\bullet$ For further details on using yfinance, see the official documentation, particularly the yf.download() API reference.

In [1]:
import yfinance as yf  # pip install yfinance if needed

raw = yf.download(
    ['TSLA', 'GS'],
    start='2013-01-01',
    end='2020-08-15',
    auto_adjust=False,
    progress=False,
    threads=False
)

prices = raw['Adj Close'][['TSLA', 'GS']]
prices.head()
Out[1]:
Ticker TSLA GS
Date
2013-01-02 2.357333 101.678612
2013-01-03 2.318000 101.122574
2013-01-04 2.293333 103.879585
2013-01-07 2.289333 103.686531
2013-01-08 2.245333 102.752068

Q2 — Simple and Log Returns¶

$\bullet$ Define the daily simple return and log return by $$ R_t=\frac{P_t}{P_{t-1}}-1,\qquad r_t=\log\!\left(\frac{P_t}{P_{t-1}}\right)=\log(1+R_t). $$

$\bullet$ Calculate both return series for each stock and store them in a new DataFrame df.

In [2]:
import numpy as np
import pandas as pd
from scipy import stats
from IPython.display import display

# Calculate daily simple and log returns
df = pd.DataFrame(index=prices.index)

for ticker in prices.columns:
    df[f'{ticker}_simple'] = prices[ticker].pct_change()
    df[f'{ticker}_log'] = np.log(prices[ticker]).diff()

df = df.dropna()
df.head()
Out[2]:
TSLA_simple TSLA_log GS_simple GS_log
Date
2013-01-03 -0.016685 -0.016826 -0.005469 -0.005484
2013-01-04 -0.010642 -0.010699 0.027264 0.026899
2013-01-07 -0.001744 -0.001746 -0.001858 -0.001860
2013-01-08 -0.019220 -0.019407 -0.009012 -0.009053
2013-01-09 -0.001187 -0.001188 0.009546 0.009500

$\bullet$ Plot histograms and boxplots of daily log returns for both stocks. Compare the empirical distributions with normal densities having the same mean and standard deviation.

In [3]:
import matplotlib.pyplot as plt
from matplotlib.ticker import PercentFormatter
from scipy import stats
from scipy.stats import norm

# Shared plotting style for the problem sheet
plot_style = {
    'font.family': 'serif',
    'font.serif': ['TeX Gyre Pagella', 'DejaVu Serif'],
    'font.size': 8,
    'axes.titlesize': 9,
    'axes.labelsize': 8,
    'xtick.labelsize': 7,
    'ytick.labelsize': 7,
    'legend.fontsize': 8,
    'pdf.fonttype': 42
}

# Daily log returns
log_data = [df[f'{t}_log'].dropna() for t in ['TSLA', 'GS']]

with plt.rc_context(plot_style):
    fig, axes = plt.subplots(2, 2, figsize=(11, 7), dpi=300)

    # (1) Histograms with fitted normal densities
    for j, (ticker, r) in enumerate(zip(['TSLA', 'GS'], log_data)):
        ax = axes[0, j]

        ax.hist(
            r, bins=60, density=True,
            color='tab:blue', alpha=0.75,
            edgecolor='white', linewidth=0.3,
            label='Observed log returns'
        )

        x = np.linspace(r.min(), r.max(), 500)
        ax.plot(
            x, norm.pdf(x, r.mean(), r.std()),
            color='tab:red', lw=1.5,
            label='Normal: same mean and SD'
        )

        ax.set_title(f'{ticker}: daily log returns')
        ax.set_xlabel('Daily log return')
        ax.set_ylabel('Density')
        ax.xaxis.set_major_formatter(PercentFormatter(1))
        ax.legend(frameon=True)
        ax.grid(alpha=0.20, color='0.85', linewidth=0.6)
        ax.set_axisbelow(True)

    # (2) Boxplots
    ax = axes[1, 0]
    ax.boxplot(
        log_data,
        tick_labels=['TSLA', 'GS'],
        patch_artist=True,
        boxprops=dict(facecolor='tab:blue', alpha=0.35),
        medianprops=dict(color='tab:orange', linewidth=1.2),
        flierprops=dict(
            marker='o', markersize=2.5,
            markerfacecolor='none',
            markeredgecolor='0.4', alpha=0.5
        )
    )

    ax.set_title('Boxplots of daily log returns')
    ax.set_ylabel('Daily log return')
    ax.yaxis.set_major_formatter(PercentFormatter(1))
    ax.grid(alpha=0.20, color='0.85', linewidth=0.6)
    ax.set_axisbelow(True)

    # (3) Normal Q-Q plot of standardized returns
    ax = axes[1, 1]

    for ticker, r, color in zip(
        ['TSLA', 'GS'], log_data, ['tab:blue', 'tab:orange']
    ):
        z = (r - r.mean()) / r.std()
        theoretical, observed = stats.probplot(
            z, dist='norm', fit=False
        )
        ax.scatter(
            theoretical, observed,
            s=7, color=color, alpha=0.65,
            label=ticker
        )

    ax.plot(
        [-4, 4], [-4, 4],
        '--', color='0.35', linewidth=1
    )

    ax.set_title('Normal Q-Q plot of standardized returns')
    ax.set_xlabel('Theoretical normal quantile')
    ax.set_ylabel('Observed standardized quantile')
    ax.legend(frameon=True)
    ax.grid(alpha=0.20, color='0.85', linewidth=0.6)
    ax.set_axisbelow(True)

    fig.subplots_adjust(
        left=0.08, right=0.98,
        bottom=0.09, top=0.94,
        wspace=0.22, hspace=0.38
    )

    plt.show()

$\bullet$ Observations. Both stocks exhibit more peaked distributions and heavier tails than the corresponding normal distributions, indicating departures from normality. The boxplots reveal numerous extreme returns, with Tesla showing greater dispersion than Goldman Sachs.

Q3 — Monte Carlo Simulation of a Trading Strategy¶

Step 1. Model Setup and Parameter Estimation

$\bullet$ Consider an initial investment of $1 million in Tesla stock on 14 Aug 2020, with a trading horizon of 100 days. The position is liquidated when the daily closing price falls below 90% or rises above 120% of the initial price, or at the end of the horizon.

$\bullet$ Assume that daily log returns are independent and normally distributed:

$$ r_t=\log\left(\frac{S_t}{S_{t-1}}\right) \sim N(\mu,\sigma^2). $$

Estimate $\mu$ and $\sigma$ using the historical daily log returns.

In [4]:
# Initial investment and trading parameters
initial_investment = 1_000_000
S0 = prices.loc['2020-08-14', 'TSLA']
n_days = 100

stop_loss = 0.90 * S0
take_profit = 1.20 * S0

# Estimate daily log-return parameters
mu_hat = df['TSLA_log'].mean()
sigma_hat = df['TSLA_log'].std()

print(f'Daily mean: {mu_hat:.6f}')
print(f'Daily standard deviation: {sigma_hat:.6f}')

Daily mean: 0.002004
Daily standard deviation: 0.034189

Step 2. Monte Carlo Simulation under Normality

$\bullet$ Simulate 100-day stock price paths using the estimated normal distribution:

$$ S_t=S_{t-1}\exp(r_t),\qquad r_t\sim N(\hat\mu,\hat\sigma^2). $$

$\bullet$ For each simulated path, liquidate the position at the first daily closing price below $0.9S_0$ or above $1.2S_0$. Otherwise, liquidate at the end of day 100.

In [5]:
# Monte Carlo simulation under normality
n_sim = 50_000
rng = np.random.default_rng(2026)

r = rng.normal(mu_hat, sigma_hat, size=(n_sim, n_days))
paths = S0 * np.exp(np.cumsum(r, axis=1))

# First stop-loss / take-profit hit, or day 100
hit = (paths < stop_loss) | (paths > take_profit)
exit_idx = np.where(hit.any(axis=1), hit.argmax(axis=1), n_days - 1)
exit_prices = paths[np.arange(n_sim), exit_idx]

# Calculate profit/loss and expected return
pnl = initial_investment * (exit_prices / S0 - 1)

loss_prob = np.mean(pnl < 0)
expected_pnl = np.mean(pnl)
expected_return = expected_pnl / initial_investment

print(f'Probability of loss: {loss_prob:.2%}')
print(f'Expected P&L (Profit / Loss): ${expected_pnl:,.2f}')
print(f'Expected simple return: {expected_return:.2%}')

# Histogram of simulated profit/loss
pnl_k = pnl / 1e3

with plt.rc_context(plot_style):
    fig, ax = plt.subplots(figsize=(7, 3.4), dpi=300)

    # Highlight regions beyond the trading thresholds
    ax.axvspan(pnl_k.min()-10, -100, color='0.85', alpha=0.25)
    ax.axvspan(200, pnl_k.max()+10, color='tab:orange', alpha=0.08)

    # Align histogram bins with trading thresholds
    bin_width = 5  # $5,000
    bin_start = -100 + bin_width * np.floor((pnl_k.min() + 100) / bin_width)
    bin_end = -100 + bin_width * np.ceil((pnl_k.max() + 100) / bin_width)
    bins = np.arange(bin_start, bin_end + bin_width, bin_width)

    ax.hist(
        pnl_k, bins=bins,
        color='tab:blue', alpha=0.88,
        edgecolor='white', linewidth=0.35
    )

    ax.axvline(
        -100, color='0.35', linestyle='--',
        linewidth=1.1, label='Stop-loss threshold (-10%)'
    )
    ax.axvline(
        200, color='tab:orange', linestyle='--',
        linewidth=1.1, label='Take-profit threshold (+20%)'
    )

    ax.set_xlabel('Profit / Loss ($ thousands)')
    ax.set_ylabel('Frequency')
    ax.set_title('Profit / Loss under Normal Log Returns')
    ax.legend(frameon=False, fontsize=8)

    ax.set_xlim(pnl_k.min()-10, pnl_k.max()+10)
    ax.set_axisbelow(True)
    ax.grid(alpha=0.20, color='0.85', linewidth=0.6)

    fig.tight_layout()
    plt.show()

Probability of loss: 48.25%
Expected P&L (Profit / Loss): $59,457.90
Expected simple return: 5.95%

Step 3. Monte Carlo Simulation under Student-t Log Returns

$\bullet$ Replace the normal distribution with a Student-t distribution with 7 degrees of freedom:

$$ r_t \sim t_7(\mathrm{loc},\mathrm{scale}). $$

$\bullet$ Estimate the location and scale parameters using scipy.stats.t.fit, fixing the degrees of freedom at 7.

$\bullet$ Repeat the Monte Carlo simulation under the same trading strategy. Calculate the probability of loss, expected P&L, and expected simple return, and plot the P&L distribution.

$\bullet$ Student-t probability density function is

$$ f(r;\nu,\ell,s) = \frac{\Gamma\left(\frac{\nu+1}{2}\right)} {s\sqrt{\nu\pi}\,\Gamma\left(\frac{\nu}{2}\right)} \left(1+\frac{(r-\ell)^2}{\nu s^2}\right)^{-\frac{\nu+1}{2}}, $$

where $\nu=7$ is the degrees of freedom, $\ell$ is the location parameter, and $s>0$ is the scale parameter.

$\bullet$ scipy.stats.t.fit estimates $\ell$ and $s$ by maximum likelihood.

In [6]:
# Fit Student-t distribution with fixed degrees of freedom
nu = 7
_, loc_t, scale_t = stats.t.fit(df['TSLA_log'], fdf=nu)

print(f'Degrees of freedom: {nu}')
print(f'Location: {loc_t:.6f}')
print(f'Scale: {scale_t:.6f}')

# Monte Carlo simulation under Student-t
rng_t = np.random.default_rng(2027)
r_t = loc_t + scale_t * rng_t.standard_t(
    nu, size=(n_sim, n_days)
)
paths_t = S0 * np.exp(np.cumsum(r_t, axis=1))

# Liquidate at the first barrier crossing or day 100
hit_t = (paths_t < stop_loss) | (paths_t > take_profit)
exit_idx_t = np.where(
    hit_t.any(axis=1), hit_t.argmax(axis=1), n_days - 1
)
exit_prices_t = paths_t[np.arange(n_sim), exit_idx_t]

# Profit/loss and expected return
pnl_t = initial_investment * (exit_prices_t / S0 - 1)

loss_prob_t = np.mean(pnl_t < 0)
expected_pnl_t = np.mean(pnl_t)
expected_return_t = expected_pnl_t / initial_investment

print(f'\nProbability of loss: {loss_prob_t:.2%}')
print(f'Expected P&L: ${expected_pnl_t:,.2f}')
print(f'Expected simple return: {expected_return_t:.2%}')

# Histogram of simulated profit/loss
pnl_t_k = pnl_t / 1e3

with plt.rc_context(plot_style):
    fig, ax = plt.subplots(figsize=(7, 3.4), dpi=300)

    ax.axvspan(pnl_t_k.min()-10, -100, color='0.85', alpha=0.25)
    ax.axvspan(200, pnl_t_k.max()+10, color='tab:orange', alpha=0.08)

    # Align bins with trading thresholds
    bin_width = 5
    bin_start = -100 + bin_width * np.floor((pnl_t_k.min()+100) / bin_width)
    bin_end = -100 + bin_width * np.ceil((pnl_t_k.max()+100) / bin_width)
    bins = np.arange(bin_start, bin_end + bin_width, bin_width)

    ax.hist(
        pnl_t_k, bins=bins,
        color='tab:blue', alpha=0.88,
        edgecolor='white', linewidth=0.35
    )

    ax.axvline(
        -100, color='0.35', linestyle='--',
        linewidth=1.1, label='Stop-loss threshold (-10%)'
    )
    ax.axvline(
        200, color='tab:orange', linestyle='--',
        linewidth=1.1, label='Take-profit threshold (+20%)'
    )

    ax.set_xlabel('Profit / Loss ($ thousands)')
    ax.set_ylabel('Frequency')
    ax.set_title('Profit / Loss under Student-t Log Returns')
    ax.legend(frameon=False, fontsize=8)

    ax.set_xlim(pnl_t_k.min()-10, pnl_t_k.max()+10)
    ax.set_axisbelow(True)
    ax.grid(alpha=0.20, color='0.85', linewidth=0.6)

    fig.tight_layout()
    plt.show()

Degrees of freedom: 7
Location: 0.001636
Scale: 0.025173

Probability of loss: 47.59%
Expected P&L: $60,180.72
Expected simple return: 6.02%

Step 4. Comparison of Normal and Student-t Models

$\bullet$ Compare the probability of loss, average P&L, and average simple return under the two models. Examine the extreme outcomes using the 1st and 99th percentiles, together with the minimum and maximum simulated P&L.

In [7]:
# Compare Monte Carlo results
comparison = pd.DataFrame({
    'Normal': {
        'Probability of loss': np.mean(pnl < 0),
        'Mean P&L ($)': np.mean(pnl),
        'Mean return': np.mean(pnl) / initial_investment,
        '1% P&L ($)': np.percentile(pnl, 1),
        '99% P&L ($)': np.percentile(pnl, 99),
        'Min P&L ($)': np.min(pnl),
        'Max P&L ($)': np.max(pnl)
    },
    'Student-t': {
        'Probability of loss': np.mean(pnl_t < 0),
        'Mean P&L ($)': np.mean(pnl_t),
        'Mean return': np.mean(pnl_t) / initial_investment,
        '1% P&L ($)': np.percentile(pnl_t, 1),
        '99% P&L ($)': np.percentile(pnl_t, 99),
        'Min P&L ($)': np.min(pnl_t),
        'Max P&L ($)': np.max(pnl_t)
    }
}).T


# Format comparison table without additional dependencies
comparison_display = comparison.copy()

# Percentage columns
for col in ['Probability of loss', 'Mean return']:
    comparison_display[col] = comparison_display[col].map(
        lambda x: f'{x:.2%}'
    )

# Monetary columns
for col in comparison.columns.difference(
    ['Probability of loss', 'Mean return']
):
    comparison_display[col] = comparison_display[col].map(
        lambda x: f'{x:,.0f}'
    )

display(comparison_display)
Probability of loss Mean P&L ($) Mean return 1% P&L ($) 99% P&L ($) Min P&L ($) Max P&L ($)
Normal 48.25% 59,458 5.95% -155,985 284,431 -199,548 353,341
Student-t 47.59% 60,181 6.02% -162,543 290,552 -313,199 612,305

$\bullet$ Observations. The Student-t model produces more extreme simulated gains and losses than the normal model, reflecting its heavier tails. The reported mean P&L is a finite-simulation average and is sensitive to extreme outcomes (since Student-t model may not have corresponding P\&L expecation).

Part B¶

Q1–Q2 — Load and Explore the Fuel Data¶

$\bullet$ Load FuelData.csv into a Pandas DataFrame using pd.read_csv().

$\bullet$ Inspect the dataset using basic Pandas operations, including its dimensions, variable types, and summary statistics.

In [8]:
# Load the fuel efficiency dataset
fuel = pd.read_csv('FuelData.csv')

# Inspect the dataset
fuel.info()
display(fuel.head())
display(fuel.describe())

<class 'pandas.DataFrame'>
RangeIndex: 1234 entries, 0 to 1233
Columns: 167 entries, Model Year to Unnamed: 166
dtypes: float64(82), int64(19), str(66)
memory usage: 1.6 MB
Model Year Mfr Name Division Carline Verify Mfr Cd Index (Model Type Index) Eng Displ # Cyl Transmission City FE (Guide) - Conventional Fuel ... 120V Charge time at 120 Volts (hours) PHEV Total Driving Range (rounded to nearest 10 miles)DISTANCE City PHEV Composite MPGe Hwy PHEV Composite MPGe Comb PHEV Composite MPGe Unnamed: 162 Unnamed: 163 Unnamed: 164 Unnamed: 165 Unnamed: 166
0 2019 Honda Acura NSX HNX 41 3.5 6 Auto(AM-S9) 21 ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
1 2019 FCA US LLC ALFA ROMEO 4C CRX 403 1.8 4 Auto(AM6) 24 ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2 2019 aston martin Aston Martin Lagonda Ltd Vantage V8 ASX 3 4.0 8 Auto(S8) 18 ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
3 2019 Volkswagen Group of Audi TT Roadster quattro VGA 51 2.0 4 Auto(AM-S7) 23 ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
4 2019 BMW BMW Z4 sDrive30i BMX 350 2.0 4 Auto(S8) 25 ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN

5 rows × 167 columns

Model Year Index (Model Type Index) Eng Displ # Cyl City FE (Guide) - Conventional Fuel Hwy FE (Guide) - Conventional Fuel Comb FE (Guide) - Conventional Fuel City Unadj FE - Conventional Fuel Hwy Unadj FE - Conventional Fuel Comb Unadj FE - Conventional Fuel ... 120V Charge time at 120 Volts (hours) PHEV Total Driving Range (rounded to nearest 10 miles)DISTANCE City PHEV Composite MPGe Hwy PHEV Composite MPGe Comb PHEV Composite MPGe Unnamed: 162 Unnamed: 163 Unnamed: 164 Unnamed: 165 Unnamed: 166
count 1234.0 1234.000000 1234.000000 1234.000000 1234.000000 1234.000000 1234.000000 1234.000000 1234.000000 1234.000000 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
mean 2019.0 222.199352 3.067666 5.535656 20.760940 27.701783 23.326580 26.536126 39.685080 31.075385 ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
std 0.0 225.366467 1.294934 1.774134 6.040716 6.035866 5.963758 8.685933 9.141559 8.707823 ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
min 2019.0 1.000000 1.000000 3.000000 9.000000 14.000000 11.000000 10.400000 18.200000 12.885000 ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
25% 2019.0 41.000000 2.000000 4.000000 17.000000 24.000000 19.000000 21.000000 33.200000 25.092475 ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
50% 2019.0 126.000000 3.000000 6.000000 20.000000 27.000000 23.000000 25.026100 38.700000 29.845800 ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
75% 2019.0 385.750000 3.600000 6.000000 23.000000 31.000000 26.000000 29.827700 44.454600 35.023200 ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
max 2019.0 991.000000 8.000000 16.000000 58.000000 59.000000 58.000000 83.559800 79.100000 80.847200 ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN

8 rows × 101 columns

Q3–Q4 — Fuel Cost and Efficiency Distributions¶

$\bullet$ Plot a histogram of annual fuel cost for conventional-fuel vehicles.

$\bullet$ Plot boxplots of combined fuel efficiency (conventional fuel), grouped by manufacturer, to compare their distributions.

In [9]:
# Q3–Q4: Fuel cost and efficiency distributions
cost = fuel['Annual Fuel1 Cost - Conventional Fuel'].dropna()

col = 'Comb FE (Guide) - Conventional Fuel'
groups = fuel.groupby('Mfr Name')[col]
order = groups.median().sort_values().index
data = [groups.get_group(name).dropna() for name in order]

with plt.rc_context(plot_style):
    fig, axes = plt.subplots(
        1, 2, figsize=(12, 5.2), dpi=300,
        gridspec_kw={'width_ratios': [1.3, 1]}
    )

    # (a) Annual fuel cost
    ax = axes[0]
    ax.hist(
        cost, bins=25,
        color='tab:blue', alpha=0.88,
        edgecolor='white', linewidth=0.35
    )
    ax.set_xlabel('Annual Fuel Cost ($)')
    ax.set_ylabel('Frequency')
    ax.set_title('(a) Annual Fuel Cost')

    # (b) Fuel efficiency by manufacturer
    ax = axes[1]
    ax.boxplot(
        data,
        orientation='horizontal',
        patch_artist=True,
        widths=0.6,
        boxprops=dict(facecolor='tab:blue', alpha=0.4),
        medianprops=dict(color='tab:orange', linewidth=1.2),
        flierprops=dict(
            marker='o', markersize=2,
            markerfacecolor='none',
            markeredgecolor='0.4', alpha=0.5
        )
    )

    ax.set_yticks(np.arange(1, len(order) + 1))
    ax.set_yticklabels(order, fontsize=6.5)
    ax.set_xlabel('Combined Fuel Efficiency (MPG)')
    ax.set_title('(b) Fuel Efficiency by Manufacturer')

    # Shared style
    for ax in axes:
        ax.set_axisbelow(True)
        ax.grid(alpha=0.20, color='0.85', linewidth=0.6)

    fig.subplots_adjust(
        left=0.065, right=0.985,
        bottom=0.12, top=0.92,
        wspace=0.40
    )
    plt.show()

Q5 — Most and Least Fuel-Efficient Cars¶

$\bullet$ Identify the cars with the highest and lowest combined fuel efficiency (MPG) in the dataset, including any ties.

In [10]:
# Identify the most and least fuel-efficient cars
col = 'Comb FE (Guide) - Conventional Fuel'
cols = ['Mfr Name', 'Carline', col]

print('Most fuel-efficient cars:')
display(fuel.loc[fuel[col] == fuel[col].max(), cols])

print('Least fuel-efficient cars:')
display(fuel.loc[fuel[col] == fuel[col].min(), cols])

Most fuel-efficient cars:
Mfr Name Carline Comb FE (Guide) - Conventional Fuel
634 Hyundai Ioniq Blue 58 
Least fuel-efficient cars:
Mfr Name Carline Comb FE (Guide) - Conventional Fuel
5 Volkswagen Group of Chiron 11
38 Volkswagen Group of Aventador Coupe 11
39 Volkswagen Group of Aventador Roadster 11

Q6–Q7 — Fuel Efficiency and Engine Characteristics¶

$\bullet$ Plot combined fuel efficiency against (a) engine displacement and (b) the number of cylinders.

In [11]:
# Fuel efficiency and engine characteristics
y = fuel['Comb FE (Guide) - Conventional Fuel']
x1 = fuel['Eng Displ']
x2 = fuel['# Cyl']

rng = np.random.default_rng(2026)

with plt.rc_context(plot_style):
    fig, axes = plt.subplots(1, 2, figsize=(10, 3.5), dpi=300)

    # (a) Engine displacement
    axes[0].scatter(
        x1, y, s=15, color='tab:blue',
        alpha=0.45, edgecolors='none'
    )
    axes[0].set_xlabel('Engine Displacement (L)')
    axes[0].set_ylabel('Combined Fuel Efficiency (MPG)')
    axes[0].set_title('(a) Engine Displacement')

    # (b) Number of cylinders
    x2_jitter = x2 + rng.uniform(-0.15, 0.15, len(x2))
    axes[1].scatter(
        x2_jitter, y, s=15, color='tab:blue',
        alpha=0.35, edgecolors='none'
    )
    axes[1].set_xticks(sorted(x2.dropna().unique()))
    axes[1].set_xlabel('Number of Cylinders')
    axes[1].set_title('(b) Number of Cylinders')

    for ax in axes:
        ax.set_axisbelow(True)
        ax.grid(alpha=0.20, color='0.85', linewidth=0.6)

    fig.tight_layout()
    plt.show()

print(f'Correlation (Displacement): {x1.corr(y):.3f}')
print(f'Correlation (Cylinders):    {x2.corr(y):.3f}')

Correlation (Displacement): -0.720
Correlation (Cylinders):    -0.688

Q8–Q9 — Simple Linear Regression¶

$\bullet$ Define the variables:

  • $x$: Engine Displacement (L).
  • $y$: Combined Fuel Efficiency (MPG).
  • $z$: Transmission Type (Automatic / Manual), treated as a categorical variable.

$\bullet$ Estimate the simple linear regression model using ordinary least squares (OLS):

$$ y_i=\beta_0+\beta_1 x_i+\varepsilon_i. $$

Report the regression summary and plot the fitted regression line.

In [12]:
import statsmodels.formula.api as smf

# Define regression variables
reg = pd.DataFrame({
    'x': fuel['Eng Displ'],
    'y': fuel['Comb FE (Guide) - Conventional Fuel'],
    'z': pd.Categorical(
        np.where(
            fuel['Transmission'].str.startswith('Manual'),
            'Manual', 'Automatic'
        ),
        categories=['Automatic', 'Manual']
    )
}).dropna()

# Simple linear regression
model1 = smf.ols('y ~ x', data=reg).fit()
print(model1.summary())

# Scatterplot with fitted regression line
x_grid = np.linspace(reg['x'].min(), reg['x'].max(), 200)
y_fit = model1.predict(pd.DataFrame({'x': x_grid}))

with plt.rc_context(plot_style):
    fig, ax = plt.subplots(figsize=(7, 3.4), dpi=300)

    ax.scatter(
        reg['x'], reg['y'],
        s=15, color='tab:blue', alpha=0.4,
        edgecolors='none', label='Observed'
    )
    ax.plot(
        x_grid, y_fit,
        color='tab:orange', linewidth=1.5,
        label='OLS fit'
    )

    ax.set_xlabel('Engine Displacement (L)')
    ax.set_ylabel('Combined Fuel Efficiency (MPG)')
    ax.set_title('Simple Linear Regression')
    ax.legend(frameon=False)

    ax.set_axisbelow(True)
    ax.grid(alpha=0.20, color='0.85', linewidth=0.6)

    fig.tight_layout()
    plt.show()

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.519
Model:                            OLS   Adj. R-squared:                  0.518
Method:                 Least Squares   F-statistic:                     1329.
Date:                Tue, 13 Oct 2026   Prob (F-statistic):          6.06e-198
Time:                        14:13:08   Log-Likelihood:                -3502.6
No. Observations:                1234   AIC:                             7009.
Df Residuals:                    1232   BIC:                             7020.
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept     33.5032      0.303    110.558      0.000      32.909      34.098
x             -3.3174      0.091    -36.449      0.000      -3.496      -3.139
==============================================================================
Omnibus:                      765.063   Durbin-Watson:                   0.953
Prob(Omnibus):                  0.000   Jarque-Bera (JB):             9521.553
Skew:                           2.688   Prob(JB):                         0.00
Kurtosis:                      15.501   Cond. No.                         9.23
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

$\bullet$ Observations. The estimated regression is

$$ \hat y = 33.503 - 3.317x, \qquad R^2 = 0.519. $$

A one-litre increase in engine displacement is associated with approximately 3.32 fewer miles per gallon. Engine displacement explains about 51.9% of the variation in combined fuel efficiency.

Q10 — Regression Diagnostics¶

$\bullet$ Regression diagnostics. The residual is the difference between observed and predicted fuel efficiency: $e_i=y_i-\hat y_i$. The orange LOWESS curves show smoothed trends in the residuals.

Assumption What to check Observation
Linearity Residuals vs $x$ should show no systematic pattern. The curved trend suggests nonlinearity.
Homoscedasticity Residual spread should remain approximately constant. The changing spread suggests heteroscedasticity.
Normality Q–Q points should follow the diagonal line. Strong upper-tail deviations suggest non-normality.
Independence Errors should be unrelated across observations. Cannot be verified from the observation index alone.

$\bullet$ Conclusion. The diagnostics suggest that the simple linear model is inadequate. A nonlinear specification or a transformation of $y$ may improve the model.

In [13]:
from scipy import stats
from statsmodels.nonparametric.smoothers_lowess import lowess

# Fitted values and residuals
fitted = model1.fittedvalues.to_numpy()
resid = model1.resid.to_numpy()
std_resid = resid / resid.std(ddof=1)

with plt.rc_context(plot_style):
    fig, axes = plt.subplots(2, 2, figsize=(10, 6.5), dpi=300)

    # (a) Linearity
    ax = axes[0, 0]
    ax.scatter(reg['x'], resid, s=12, alpha=0.4,
               color='tab:blue', edgecolors='none')
    smooth = lowess(resid, reg['x'], frac=0.3)
    ax.plot(smooth[:, 0], smooth[:, 1],
            color='tab:orange', lw=1.5)
    ax.axhline(0, color='0.4', ls='--', lw=0.9)
    ax.set_title('(a) Linearity')
    ax.set_xlabel('Engine Displacement (L)')
    ax.set_ylabel('Residuals')

    # (b) Homoscedasticity
    ax = axes[0, 1]
    scale_resid = np.sqrt(np.abs(std_resid))
    ax.scatter(fitted, scale_resid, s=12, alpha=0.4,
               color='tab:blue', edgecolors='none')
    smooth = lowess(scale_resid, fitted, frac=0.3)
    ax.plot(smooth[:, 0], smooth[:, 1],
            color='tab:orange', lw=1.5)
    ax.set_title('(b) Homoscedasticity')
    ax.set_xlabel('Fitted Values')
    ax.set_ylabel(r'$\sqrt{|\mathrm{Standardized\ Residuals}|}$')

    # (c) Normality
    ax = axes[1, 0]
    theoretical, observed = stats.probplot(
        std_resid, dist='norm', fit=False
    )
    ax.scatter(theoretical, observed, s=10, alpha=0.5,
               color='tab:blue', edgecolors='none')
    limits = [
        min(min(theoretical), min(observed)),
        max(max(theoretical), max(observed))
    ]
    ax.plot(limits, limits, '--', color='0.4', lw=0.9)
    ax.set_title('(c) Normal Q-Q Plot')
    ax.set_xlabel('Theoretical Quantiles')
    ax.set_ylabel('Observed Quantiles')

    # (d) Independence
    ax = axes[1, 1]
    ax.scatter(np.arange(len(resid)), resid, s=10,
               alpha=0.4, color='tab:blue', edgecolors='none')
    ax.axhline(0, color='0.4', ls='--', lw=0.9)
    ax.set_title('(d) Independence')
    ax.set_xlabel('Observation Index')
    ax.set_ylabel('Residuals')

    for ax in axes.flat:
        ax.set_axisbelow(True)
        ax.grid(alpha=0.20, color='0.85', linewidth=0.6)

    fig.tight_layout()
    plt.show()

Q11 — Log Transformation and Regression Diagnostics¶

$\bullet$ Repeat the regression using the logarithm of fuel efficiency as the dependent variable:

$$ \log(y_i)=\beta_0+\beta_1x_i+\varepsilon_i. $$

$\bullet$ Report the OLS summary and repeat the four diagnostic plots from Q10. Examine whether the log transformation improves the model assumptions.

In [14]:
# Regression with log-transformed fuel efficiency
reg['log_y'] = np.log(reg['y'])
model_log = smf.ols('log_y ~ x', data=reg).fit()

print(model_log.summary())

# Fitted values and residuals
fitted = model_log.fittedvalues.to_numpy()
resid = model_log.resid.to_numpy()
std_resid = resid / resid.std(ddof=1)

with plt.rc_context(plot_style):
    fig, axes = plt.subplots(2, 2, figsize=(10, 6.5), dpi=300)

    # (a) Linearity
    ax = axes[0, 0]
    ax.scatter(reg['x'], resid, s=12, alpha=0.4,
               color='tab:blue', edgecolors='none')
    smooth = lowess(resid, reg['x'], frac=0.3)
    ax.plot(smooth[:, 0], smooth[:, 1],
            color='tab:orange', lw=1.5)
    ax.axhline(0, color='0.4', ls='--', lw=0.9)
    ax.set_title('(a) Linearity')
    ax.set_xlabel('Engine Displacement (L)')
    ax.set_ylabel('Residuals')

    # (b) Homoscedasticity
    ax = axes[0, 1]
    scale_resid = np.sqrt(np.abs(std_resid))
    ax.scatter(fitted, scale_resid, s=12, alpha=0.4,
               color='tab:blue', edgecolors='none')
    smooth = lowess(scale_resid, fitted, frac=0.3)
    ax.plot(smooth[:, 0], smooth[:, 1],
            color='tab:orange', lw=1.5)
    ax.set_title('(b) Homoscedasticity')
    ax.set_xlabel('Fitted Log Fuel Efficiency')
    ax.set_ylabel(r'$\sqrt{|\mathrm{Standardized\ Residuals}|}$')

    # (c) Normality
    ax = axes[1, 0]
    theoretical, observed = stats.probplot(
        std_resid, dist='norm', fit=False
    )
    ax.scatter(theoretical, observed, s=10, alpha=0.5,
               color='tab:blue', edgecolors='none')
    limits = [
        min(min(theoretical), min(observed)),
        max(max(theoretical), max(observed))
    ]
    ax.plot(limits, limits, '--', color='0.4', lw=0.9)
    ax.set_title('(c) Normal Q-Q Plot')
    ax.set_xlabel('Theoretical Quantiles')
    ax.set_ylabel('Observed Quantiles')

    # (d) Independence
    ax = axes[1, 1]
    ax.scatter(np.arange(len(resid)), resid, s=10,
               alpha=0.4, color='tab:blue', edgecolors='none')
    ax.axhline(0, color='0.4', ls='--', lw=0.9)
    ax.set_title('(d) Independence')
    ax.set_xlabel('Observation Index')
    ax.set_ylabel('Residuals')

    for ax in axes.flat:
        ax.set_axisbelow(True)
        ax.grid(alpha=0.20, color='0.85', linewidth=0.6)

    fig.tight_layout()
    plt.show()

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                  log_y   R-squared:                       0.638
Model:                            OLS   Adj. R-squared:                  0.638
Method:                 Least Squares   F-statistic:                     2170.
Date:                Tue, 13 Oct 2026   Prob (F-statistic):          4.93e-274
Time:                        14:13:09   Log-Likelihood:                 648.99
No. Observations:                1234   AIC:                            -1294.
Df Residuals:                    1232   BIC:                            -1284.
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept      3.5702      0.010    340.666      0.000       3.550       3.591
x             -0.1466      0.003    -46.587      0.000      -0.153      -0.140
==============================================================================
Omnibus:                      317.428   Durbin-Watson:                   0.959
Prob(Omnibus):                  0.000   Jarque-Bera (JB):             1151.059
Skew:                           1.213   Prob(JB):                    1.12e-250
Kurtosis:                       7.062   Cond. No.                         9.23
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

$\bullet$ Observations. The log transformation improves the approximate normality of residuals, but the diagnostics still indicate:

  • Nonlinearity: The smoothed residual trend remains curved.
  • Heteroscedasticity: Residual dispersion is not constant.
  • Non-normality: The Q–Q plot still shows substantial upper-tail deviations.
  • Independence: Cannot be established from the observation order alone.

$\bullet$ Conclusion. The log transformation alone does not resolve the model's limitations. Adding a quadratic term in engine displacement may better capture the nonlinear relationship.

Q12 — Model Extensions and Comparison¶

$\bullet$ Examine whether adding transmission type ($z$) and a quadratic term in engine displacement ($x$) improves the regression model:

$$ \begin{aligned} M_1:\quad & y_i=\beta_0+\beta_1x_i+\varepsilon_i,\\ M_2:\quad & y_i=\beta_0+\beta_1x_i+\beta_2z_i+\varepsilon_i,\\ M_3:\quad & y_i=\beta_0+\beta_1x_i+\beta_2x_i^2+\varepsilon_i,\\ M_4:\quad & y_i=\beta_0+\beta_1x_i+\beta_2x_i^2+\beta_3z_i+\varepsilon_i. \end{aligned} $$

Here $z_i=1$ for Manual and $0$ for Automatic transmission.

$\bullet$ Compare model fit using $R^2$, adjusted $R^2$, and AIC. Use F-tests to assess whether the additional predictors significantly improve the model.

$\bullet$ F-test. To test whether additional predictors improve the model, the F-statistic is

$$ F = \frac{(RSS_R-RSS_U)/(k_U-k_R)}{RSS_U/(n-k_U)} \overset{H_0}{\sim} F_{k_U-k_R,\;n-k_U}, $$

where $RSS_R$ and $RSS_U$ are the residual sums of squares of Model $R$ and Model $U$, $k_R$ and $k_U$ are the numbers of estimated parameters (including the intercept), and $n$ is the sample size.

In [15]:
# Fit the four regression models
model2 = smf.ols('y ~ x + C(z)', data=reg).fit()
model3 = smf.ols('y ~ x + I(x**2)', data=reg).fit()
model4 = smf.ols('y ~ x + I(x**2) + C(z)', data=reg).fit()

models = {
    'M1: Linear': model1,
    'M2: + Transmission': model2,
    'M3: + Quadratic': model3,
    'M4: + Both': model4
}


model_comparison = pd.DataFrame({
    name: {
        'R²': m.rsquared,
        'Adj. R²': m.rsquared_adj,
        'AIC': m.aic
    }
    for name, m in models.items()
}).T

# Model extensions for log fuel efficiency
log_models = {
    'M1: Linear': model_log,
    'M2: + Transmission': smf.ols('log_y ~ x + C(z)', data=reg).fit(),
    'M3: + Quadratic': smf.ols('log_y ~ x + I(x**2)', data=reg).fit(),
    'M4: + Both': smf.ols('log_y ~ x + I(x**2) + C(z)', data=reg).fit()
}

log_comparison = pd.DataFrame({
    name: {
        'R²': m.rsquared,
        'Adj. R²': m.rsquared_adj,
        'AIC': m.aic
    }
    for name, m in log_models.items()
}).T

print('Panel A: Fuel Efficiency')
display(model_comparison.round(3))

print('Panel B: Log Fuel Efficiency')
display(log_comparison.round(3))


from scipy.stats import f

# Nested F-test using residual sums of squares
def nested_f_test(restricted, unrestricted):
    rss_r = restricted.ssr
    rss_u = unrestricted.ssr

    k_r = int(restricted.df_model + 1)
    k_u = int(unrestricted.df_model + 1)
    n = int(unrestricted.nobs)

    df1 = k_u - k_r
    df2 = n - k_u

    F_stat = ((rss_r - rss_u) / df1) / (rss_u / df2)
    p_value = f.sf(F_stat, df1, df2)

    return F_stat, p_value

# Compare nested models within each panel
pairs = [
    ('M1 vs M2', 'M1: Linear', 'M2: + Transmission'),
    ('M1 vs M3', 'M1: Linear', 'M3: + Quadratic'),
    ('M2 vs M4', 'M2: + Transmission', 'M4: + Both'),
    ('M3 vs M4', 'M3: + Quadratic', 'M4: + Both')
]

for panel, model_dict in [
    ('Panel A: Fuel Efficiency', models),
    ('Panel B: Log Fuel Efficiency', log_models)
]:
    rows = []
    for label, r, u in pairs:
        F_stat, p_value = nested_f_test(
            model_dict[r], model_dict[u]
        )
        rows.append({
            'Comparison': label,
            'F-statistic': F_stat,
            'p-value': p_value
        })

    table = pd.DataFrame(rows).set_index('Comparison')
    table['F-statistic'] = table['F-statistic'].round(3)
    table['p-value'] = table['p-value'].map(lambda p: f'{p:.3g}')

    print(panel)
    display(table)

Panel A: Fuel Efficiency
R² Adj. R² AIC
M1: Linear 0.519 0.518 7009.278
M2: + Transmission 0.520 0.520 7007.569
M3: + Quadratic 0.571 0.570 6870.341
M4: + Both 0.575 0.574 6860.633 
Panel B: Log Fuel Efficiency
R² Adj. R² AIC
M1: Linear 0.638 0.638 -1293.984
M2: + Transmission 0.639 0.638 -1295.593
M3: + Quadratic 0.667 0.666 -1395.258
M4: + Both 0.670 0.669 -1403.219 
Panel A: Fuel Efficiency
F-statistic p-value
Comparison
M1 vs M2 3.706 0.0544
M1 vs M3 148.939 2.03e-32
M2 vs M4 157.783 3.89e-34
M3 vs M4 11.725 0.000637 
Panel B: Log Fuel Efficiency
F-statistic p-value
Comparison
M1 vs M2 3.606 0.0578
M1 vs M3 107.458 3.38e-24
M2 vs M4 114.272 1.45e-25
M3 vs M4 9.969 0.00163

Q13 — Cars with Extreme Residual Fuel Efficiency¶

$\bullet$ Using Model $M_4$, identify the cars with the largest positive and negative residuals after accounting for engine displacement and transmission type:

$$ e_i=y_i-\hat y_i. $$

A positive residual indicates higher fuel efficiency than predicted, while a negative residual indicates lower fuel efficiency than predicted.

In [16]:
# Calculate residuals from the full model (M4)
results = fuel.loc[reg.index, ['Mfr Name', 'Carline']].copy()

results['Engine (L)'] = reg['x']
results['Transmission'] = reg['z']
results['Observed MPG'] = reg['y']
results['Predicted MPG'] = model4.predict(reg)
results['Residual (MPG)'] = (
    results['Observed MPG'] - results['Predicted MPG']
)

# Identify the largest positive and negative residuals
extremes = results.loc[[
    results['Residual (MPG)'].idxmax(),
    results['Residual (MPG)'].idxmin()
]].copy()

extremes.index = ['Largest positive', 'Largest negative']

display(extremes.round(2))
Mfr Name Carline Engine (L) Transmission Observed MPG Predicted MPG Residual (MPG)
Largest positive Hyundai Ioniq Blue 1.6 Automatic 58 30.01 27.99
Largest negative Volkswagen Group of Chiron 8.0 Automatic 11 19.06 -8.06