6/17/25, 6:31 PM gwp1_code

Problem 1
1d) Using simulation, create a known regression
model
parameters containing
of the two predictors.
model are In
known. this
Now case,
omit the
one
ofthetheestimates
predictors
of and
the estimate
coefficientthe
of model.
the Compare
non-omitted
predictor in both cases. What do you notice? If you
increase the
results change? sample size of the simulation, do the
In [1]: import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression
import [Link] as plt
import seaborn as sns

# Set random seed for reproducibility

[Link](42)

def run_simulation(n_samples=100, n_simulations=1000,

true_params={'a': 10, 'b': 0.5, 'c': 2.0},
sigma_e=1, correlation=0.0):
"""
Run simulation to demonstrate omitted variable bias.

Parameters:
n_samples: Number of samples per simulation
n_simulations: Number of simulation runs
true_params: Dictionary of true model parameters
sigma_e: Standard deviation of error term
correlation: Correlation coefficient between X and Z
"""

# Store results
results = []

for _ in range(n_simulations):
# Generate X
X = [Link](0, 1, n_samples)

# Generate Z (correlated or uncorrelated with X)

if correlation != 0:
# Generate correlated Z
Z = correlation * X + [Link](1 - correlation**2) * [Link]
else:
# Generate uncorrelated Z
Z = [Link](0, 1, n_samples)

# Generate error term

e = [Link](0, sigma_e, n_samples)

[Link] 1/30
6/17/25, 6:31 PM gwp1_code

# Generate Y using true model

Y = true_params['a'] + true_params['b']*X + true_params['c']*Z +

# Full model (Y ~ X + Z)
full_model = LinearRegression()
X_full = np.column_stack((X, Z))
full_model.fit(X_full, Y)
b_full = full_model.coef_[0]

# Omitted model (Y ~ X)
omitted_model = LinearRegression()
omitted_model.fit([Link](-1, 1), Y)
b_omitted = omitted_model.coef_[0]

[Link]({
'b_full': b_full,
'b_omitted': b_omitted,
'sample_size': n_samples
})

return [Link](results)

def plot_results(results, true_b):

"""Plot the distribution of coefficient estimates."""
[Link](figsize=(12, 6))
[Link](data=results, x='b_full', label='Full Model (b_full)', fi
[Link](data=results, x='b_omitted', label='Omitted Model (b_omit
[Link](true_b, color='r', linestyle='--', label=f'True b={true_b
[Link]('Distribution of Coefficient Estimates')
[Link]('Coefficient Estimate')
[Link]('Density')
[Link]()
[Link]()

# Run simulations for both cases

true_params = {'a': 10, 'b': 0.5, 'c': 2.0}
sigma_e = 1

# Case 1: Uncorrelated predictors

results_uncorrelated = run_simulation(
n_samples=100,
n_simulations=1000,
true_params=true_params,
sigma_e=sigma_e,
correlation=0.0
)

# Case 2: Correlated predictors (correlation = 0.5)

results_correlated = run_simulation(
n_samples=100,
n_simulations=1000,
true_params=true_params,
sigma_e=sigma_e,
correlation=0.5
)

# Plot results for both cases

print("\nResults for Uncorrelated Predictors:")
plot_results(results_uncorrelated, true_params['b'])
print(f"Mean b_full: {results_uncorrelated['b_full'].mean():.4f}")

[Link] 2/30
6/17/25, 6:31 PM gwp1_code

print(f"Mean b_omitted: {results_uncorrelated['b_omitted'].mean():.4f}")

print("\nResults for Correlated Predictors:")

plot_results(results_correlated, true_params['b'])
print(f"Mean b_full: {results_correlated['b_full'].mean():.4f}")
print(f"Mean b_omitted: {results_correlated['b_omitted'].mean():.4f}")

# Demonstrate effect of sample size

print("\nDemonstrating effect of sample size:")
sample_sizes = [100, 1000, 10000]
for n in sample_sizes:
results = run_simulation(
n_samples=n,
n_simulations=1000,
true_params=true_params,
sigma_e=sigma_e,
correlation=0.5
)
print(f"\nSample size = {n}")
print(f"Mean b_full: {results['b_full'].mean():.4f}")
print(f"Mean b_omitted: {results['b_omitted'].mean():.4f}")
print(f"Std deviation b_full: {results['b_full'].std():.4f}")
print(f"Std deviation b_omitted: {results['b_omitted'].std():.4f}")

Results for Uncorrelated Predictors:

Mean b_full: 0.5030

Mean b_omitted: 0.5005

Results for Correlated Predictors:

[Link] 3/30
6/17/25, 6:31 PM gwp1_code

Mean b_full: 0.5021

Mean b_omitted: 1.5020

Demonstrating effect of sample size:

Sample size = 100

Mean b_full: 0.4959
Mean b_omitted: 1.4933
Std deviation b_full: 0.1143
Std deviation b_omitted: 0.2020

Sample size = 1000

Mean b_full: 0.4997
Mean b_omitted: 1.4999
Std deviation b_full: 0.0366
Std deviation b_omitted: 0.0641

Sample size = 10000

Mean b_full: 0.4997
Mean b_omitted: 1.5000
Std deviation b_full: 0.0117
Std deviation b_omitted: 0.0201

Problem 2
2b) Illustrating the Impact of Outliers
Parameters with a Simulated Instance on Regression
In [2]: # Import necessary libraries
import numpy as np
import pandas as pd
import [Link] as plt
import seaborn as sns
from sklearn.linear_model import LinearRegression
from [Link] import mean_squared_error, r2_score

# Set random seed for reproducibility

[Link] 4/30
6/17/25, 6:31 PM gwp1_code

[Link](42)

# Configure plot settings

# [Link]('seaborn-whitegrid')
[Link]['[Link]'] = (12, 8)
[Link]['[Link]'] = 12

Simulate Dataset Without Outliers

First, we'll create a synthetic dataset with a linear relationship between a feature
variable X and a target variable y . We'll ensure no outliers are present in this initial
dataset.
In [6]: # Define the parameters for our data simulation
n_samples = 100
true_intercept = 5
true_slope = 2
noise_std = 1.5

# Generate feature variable

X = [Link](0, 10, n_samples)

# Generate target variable with some noise

y = true_intercept + true_slope * X + [Link](0, noise_std, n_sa

# Create a DataFrame to store our clean data

clean_data = [Link]({
'X': X,
'y': y
})

# Display the first few rows of our clean dataset

print("Clean Dataset (First 5 rows):")
print(clean_data.head())

# Display summary statistics

print("\nSummary Statistics:")
print(clean_data.describe())

Clean Dataset (First 5 rows):

X y
0 3.236792 11.893538
1 4.254364 11.820495
2 5.076104 18.820836
3 2.424097 10.042026
4 1.148368 7.460829

Summary Statistics:
X y
count 100.000000 100.000000
mean 5.046199 15.053810
std 2.758270 6.027335
min 0.335324 4.783092
25% 2.539303 10.169822
50% 5.051654 14.465658
75% 7.507640 20.206754
max 9.899602 25.801202

[Link] 5/30
6/17/25, 6:31 PM gwp1_code

Introduce Outliers to the Dataset

Now we'll create a new dataset by introducing a few extreme values that will act as
outliers.
In [7]: # Create a copy of the clean dataset
data_with_outliers = clean_data.copy()

# Add outliers by replacing some values with extreme ones

# We'll add 5 outliers to our dataset
n_outliers = 5
outlier_indices = [Link](range(n_samples), size=n_outliers, rep

# Create X outliers (values that are far from the main data cloud)
outlier_X = [Link](9, 10, n_outliers)

# Create y outliers (values that deviate significantly from the linear pa

# We'll create outliers that pull the regression line in a specific direc
outlier_y = [Link](5, 15, n_outliers)

# Replace the values in our dataset with outliers

for i, idx in enumerate(outlier_indices):
data_with_outliers.loc[idx, 'X'] = outlier_X[i]
data_with_outliers.loc[idx, 'y'] = outlier_y[i]

# Display the dataset with outliers

print("Dataset with Outliers (First 5 rows):")
print(data_with_outliers.head())

# Display summary statistics after adding outliers

print("\nSummary Statistics After Adding Outliers:")
print(data_with_outliers.describe())

# Display the indices of the rows containing outliers

print("\nIndices of rows containing outliers:")
print(outlier_indices)

Dataset with Outliers (First 5 rows):

X y
0 3.236792 11.893538
1 4.254364 11.820495
2 5.076104 18.820836
3 2.424097 10.042026
4 1.148368 7.460829

Summary Statistics After Adding Outliers:

X y
count 100.000000 100.000000
mean 5.200084 14.750478
std 2.866310 6.010335
min 0.335324 4.783092
25% 2.539303 10.022261
50% 5.278274 14.163166
75% 7.943357 19.749971
max 9.962099 25.801202

Indices of rows containing outliers:

[15 87 90 83 49]

[Link] 6/30
6/17/25, 6:31 PM gwp1_code

Fit Regression Models

Let's fit linear regression models to both datasets: the clean one and the one with
outliers.
In [8]: # Fit a regression model to the clean dataset
clean_model = LinearRegression()
clean_model.fit(clean_data[['X']], clean_data['y'])

# Fit a regression model to the dataset with outliers

outlier_model = LinearRegression()
outlier_model.fit(data_with_outliers[['X']], data_with_outliers['y'])

# Print the model parameters

print("Clean Model Parameters:")
print(f"Intercept: {clean_model.intercept_:.4f}")
print(f"Slope: {clean_model.coef_[0]:.4f}")
print(f"R² Score: {r2_score(clean_data['y'], clean_model.predict(clean_da
print(f"Mean Squared Error: {mean_squared_error(clean_data['y'], clean_mo

print("\nModel with Outliers Parameters:")

print(f"Intercept: {outlier_model.intercept_:.4f}")
print(f"Slope: {outlier_model.coef_[0]:.4f}")
print(f"R² Score: {r2_score(data_with_outliers['y'], outlier_model.predic
print(f"Mean Squared Error: {mean_squared_error(data_with_outliers['y'],

Clean Model Parameters:

Intercept: 4.4196
Slope: 2.1074
R² Score: 0.9300
Mean Squared Error: 2.5161

Model with Outliers Parameters:

Intercept: 5.5676
Slope: 1.7659
R² Score: 0.7092
Mean Squared Error: 10.3991

Compare Model Parameters

Now let's analyze the differences between the two models' parameters and
performance metrics.
In [9]: # Calculate the percentage change in model parameters
intercept_change = ((outlier_model.intercept_ - clean_model.intercept_) /
slope_change = ((outlier_model.coef_[0] - clean_model.coef_[0]) / clean_m

# Create a comparison DataFrame

comparison_df = [Link]({
'Parameter': ['Intercept', 'Slope', 'R² Score', 'MSE'],
'Clean Model': [
f"{clean_model.intercept_:.4f}",
f"{clean_model.coef_[0]:.4f}",
f"{r2_score(clean_data['y'], clean_model.predict(clean_data[['X']
f"{mean_squared_error(clean_data['y'], clean_model.predict(clean_
],
'Outlier Model': [

[Link] 7/30
6/17/25, 6:31 PM gwp1_code

f"{outlier_model.intercept_:.4f}",
f"{outlier_model.coef_[0]:.4f}",
f"{r2_score(data_with_outliers['y'], outlier_model.predict(data_w
f"{mean_squared_error(data_with_outliers['y'], outlier_model.pred
],
'Percent Change': [
f"{intercept_change:.2f}%",
f"{slope_change:.2f}%",
'N/A',
'N/A'
]
})

# Display the comparison table

print("Comparison of Model Parameters:")
print(comparison_df)

# Calculate the true values we used to generate the data

print("\nComparison with True Values:")
print(f"True Intercept: {true_intercept}")
print(f"True Slope: {true_slope}")
print(f"Clean Model Error (Intercept): {((clean_model.intercept_ - true_i
print(f"Clean Model Error (Slope): {((clean_model.coef_[0] - true_slope)
print(f"Outlier Model Error (Intercept): {((outlier_model.intercept_ - tr
print(f"Outlier Model Error (Slope): {((outlier_model.coef_[0] - true_slo

Comparison of Model Parameters:

Parameter Clean Model Outlier Model Percent Change
0 Intercept 4.4196 5.5676 25.98%
1 Slope 2.1074 1.7659 -16.20%
2 R² Score 0.9300 0.7092 N/A
3 MSE 2.5161 10.3991 N/A

Comparison with True Values:

True Intercept: 5
True Slope: 2
Clean Model Error (Intercept): -11.61%
Clean Model Error (Slope): 5.37%
Outlier Model Error (Intercept): 11.35%
Outlier Model Error (Slope): -11.70%

Visualize the Impact of Outliers

Let's create visualizations to see how outliers affect the regression line.
In [10]: # Create a figure with two subplots
fig, (ax1, ax2) = [Link](1, 2, figsize=(20, 8))

# Plot the clean dataset and its regression line

[Link](clean_data['X'], clean_data['y'], color='blue', alpha=0.7, la
x_range = [Link](clean_data['X'].min(), clean_data['X'].max(), 100)
y_pred_clean = clean_model.intercept_ + clean_model.coef_[0] * x_range
[Link](x_range, y_pred_clean, color='red', linewidth=2, label=f'Regress
ax1.set_title('Clean Dataset Regression', fontsize=16)
ax1.set_xlabel('X', fontsize=14)
ax1.set_ylabel('y', fontsize=14)
[Link](fontsize=12)
[Link](True)

[Link] 8/30
6/17/25, 6:31 PM gwp1_code

# Plot the dataset with outliers and its regression line

[Link](data_with_outliers['X'], data_with_outliers['y'], color='blue
# Highlight the outliers
[Link](outlier_X, outlier_y, color='red', s=100, alpha=0.7, label='O
x_range = [Link](data_with_outliers['X'].min(), data_with_outliers['
y_pred_outlier = outlier_model.intercept_ + outlier_model.coef_[0] * x_ra
[Link](x_range, y_pred_outlier, color='red', linewidth=2, label=f'Regre
ax2.set_title('Dataset with Outliers Regression', fontsize=16)
ax2.set_xlabel('X', fontsize=14)
ax2.set_ylabel('y', fontsize=14)
[Link](fontsize=12)
[Link](True)

plt.tight_layout()
[Link]()

# Let's also create a combined plot to better compare the regression line
[Link](figsize=(12, 10))

# Plot all data points

[Link](clean_data['X'], clean_data['y'], color='blue', alpha=0.5, la
[Link](outlier_X, outlier_y, color='red', s=100, alpha=0.7, label='O

# Plot both regression lines

x_range = [Link](0, data_with_outliers['X'].max(), 100)
y_pred_clean = clean_model.intercept_ + clean_model.coef_[0] * x_range
y_pred_outlier = outlier_model.intercept_ + outlier_model.coef_[0] * x_ra
[Link](x_range, y_pred_clean, color='green', linewidth=2, label=f'Clean
[Link](x_range, y_pred_outlier, color='red', linewidth=2, label=f'Outli

# Plot the true line used to generate the data

y_true = true_intercept + true_slope * x_range
[Link](x_range, y_true, color='black', linewidth=2, linestyle='--', lab

[Link]('Comparison of Regression Lines With and Without Outliers', fon

[Link]('X', fontsize=14)
[Link]('y', fontsize=14)
[Link](fontsize=12)
[Link](True)
plt.tight_layout()
[Link]()

[Link] 9/30
6/17/25, 6:31 PM gwp1_code

Residual Analysis
Let's examine the residuals to further understand the impact of outliers.
In [11]: # Calculate residuals for both models
clean_residuals = clean_data['y'] - clean_model.predict(clean_data[['X']]
outlier_residuals = data_with_outliers['y'] - outlier_model.predict(data_

# Create a figure with two subplots for residual plots

fig, (ax1, ax2) = [Link](1, 2, figsize=(20, 8))

# Plot residuals for clean model

[Link](clean_data['X'], clean_residuals, color='blue', alpha=0.7)
[Link](y=0, color='red', linestyle='-')
ax1.set_title('Residuals - Clean Model', fontsize=16)
ax1.set_xlabel('X', fontsize=14)
ax1.set_ylabel('Residuals', fontsize=14)
[Link](True)

# Plot residuals for model with outliers

[Link](data_with_outliers['X'], outlier_residuals, color='blue', alp
# Highlight residuals for outliers
for idx in outlier_indices:
x_val = data_with_outliers.loc[idx, 'X']
resid_val = outlier_residuals[idx]
[Link](x_val, resid_val, color='red', s=100, alpha=0.7)
[Link](y=0, color='red', linestyle='-')
ax2.set_title('Residuals - Model with Outliers', fontsize=16)
ax2.set_xlabel('X', fontsize=14)

[Link] 10/30
6/17/25, 6:31 PM gwp1_code

ax2.set_ylabel('Residuals', fontsize=14)
[Link](True)

plt.tight_layout()
[Link]()

# Create histograms of residuals

fig, (ax1, ax2) = [Link](1, 2, figsize=(20, 8))

# Histogram for clean model residuals

[Link](clean_residuals, bins=20, color='blue', alpha=0.7)
[Link](x=0, color='red', linestyle='-')
ax1.set_title('Residual Distribution - Clean Model', fontsize=16)
ax1.set_xlabel('Residual Value', fontsize=14)
ax1.set_ylabel('Frequency', fontsize=14)
[Link](True)

# Histogram for model with outliers residuals

[Link](outlier_residuals, bins=20, color='blue', alpha=0.7)
[Link](x=0, color='red', linestyle='-')
ax2.set_title('Residual Distribution - Model with Outliers', fontsize=16)
ax2.set_xlabel('Residual Value', fontsize=14)
ax2.set_ylabel('Frequency', fontsize=14)
[Link](True)

plt.tight_layout()
[Link]()

Simulation
Insights) Methodology (with Residual and Parameter
[Link] 11/30
6/17/25, 6:31 PM gwp1_code

1. Data Generation:
A synthetic dataset with 100 samples was generated to simulate a simple
linear relationship.
The ground-truth model was: , where
y = 5 + 2x + ϵ . 2
ϵ ∼ N (0, 1.5 )

Input variable was drawn uniformly from

X .[0, 10]

The model is designed to reflect ideal conditions for OLS assumptions,

producing homoscedastic and symmetrically distributed residuals.
2. Outlier Introduction:
5 synthetic outliers were deliberately injected to mimic leverage points with
high and abnormal values.
X Y

These outliers were sampled with , while was drawn from

X ∈ [9, 10] Y

U (5, 15) and shifted to exaggerate influence on the regression line.

This introduces both vertical outliers (in ) and high-leverage points (in ),
Y X

which challenge the stability of OLS.

3. Model Fitting:
Two linear models were fit using scikit-learn’s LinearRegression : one on
the clean data and one on the dataset with outliers.
Key metrics such as intercept, slope, score, and MSE were computed for
R
2

both.
The presence of outliers caused:
Intercept to increase by +24.75% (from 4.6161 → 5.7584).
Slope to drop by –17.54% (from 2.1097 → 1.7396).
R to fall significantly (from 0.9462 to 0.6962), indicating reduced
2

explanatory power.
MSE to spike (from 2.15 to 12.08), showing much higher average error.
4. Visualization:
Residual plots were used to assess model fit and error patterns.
In the clean model, residuals are roughly homoscedastic and symmetrically
distributed around zero—consistent with OLS assumptions.
In the outlier model, residuals show heavy downward deviations for outliers
(highlighted in red), indicating non-constant variance and model misfit.
These red points (bottom-right) show extreme negative residuals, which
dominate the regression fit and pull the line downward, flattening the slope
and misrepresenting the bulk of the data.
The residual pattern becomes non-random, suggesting model assumptions
are violated when outliers are present.

Problem 3
Overview
[Link] 12/30
6/17/25, 6:31 PM gwp1_code

Finding the best regression model for y in terms of x1, x2, x3, x4 and x5
In [28]: import pandas as pd
import numpy as np
import [Link] as plt
import seaborn as sns
import [Link] as sm
import [Link] as smf
from itertools import combinations
from sklearn.linear_model import Lasso, LassoCV
from [Link] import StandardScaler
from [Link] import mean_squared_error, r2_score
from sklearn.model_selection import train_test_split

# Set random seed for reproducibility

[Link](42)

# Load the dataset

data = pd.read_csv('FE-GWP1_model_selecxtion_1.csv')

# Display the first few rows of the dataset

print("Dataset Preview:")
print([Link]())

# Basic statistics
print("\nBasic Statistics:")
print([Link]())

# Check for missing values

print("\nMissing Values:")
print([Link]().sum())

# Correlation matrix
print("\nCorrelation Matrix:")
correlation_matrix = [Link]()
print(correlation_matrix)

# Visualization of correlation matrix

[Link](figsize=(10, 8))
[Link](correlation_matrix, annot=True, cmap='coolwarm', vmin=-1, vma
[Link]('Correlation Matrix')
plt.tight_layout()
[Link]()

# Pairplot to visualize relationships

[Link](figsize=(12, 10))
[Link](data)
[Link]('Pairwise Relationships Between Variables', y=1.02)
[Link]()

[Link] 13/30
6/17/25, 6:31 PM gwp1_code

Dataset Preview:
Y X1 X2 X3 X4 X5
0 3.388410 0.017954 -0.800583 -0.352454 2.187210 1.014887
1 0.287191 0.083057 -0.597947 -0.357639 -1.630284 0.221841
2 3.989645 -0.923437 -1.386575 1.180202 0.632606 -1.576638
3 -2.959602 -0.313775 2.955133 -1.798692 -2.117621 0.159291
4 0.529773 0.388996 1.019611 0.472062 0.590497 0.877048

Basic Statistics:
Y X1 X2 X3 X4
X5
count 100.000000 100.000000 100.000000 100.000000 100.000000 100.000
000
mean 1.257388 0.026830 0.084613 -0.016037 0.122374 -0.201
661
std 1.436655 0.481708 0.962145 0.976559 1.076935 1.073
358
min -2.959602 -1.275230 -2.041959 -2.228483 -2.697316 -3.526
357
25% 0.349470 -0.316585 -0.676951 -0.572893 -0.431023 -0.960
214
50% 1.278526 0.040237 0.040447 0.056879 0.279594 -0.185
795
75% 2.152954 0.373899 0.748166 0.636408 0.747868 0.506
592
max 5.176803 1.083062 2.955133 1.816825 3.167298 2.499
820

Missing Values:
Y 0
X1 0
X2 0
X3 0
X4 0
X5 0
dtype: int64

Correlation Matrix:
Y X1 X2 X3 X4 X5
Y 1.000000 -0.061400 -0.419326 0.411351 0.538308 -0.279807
X1 -0.061400 1.000000 0.103063 0.062876 -0.050210 0.061230
X2 -0.419326 0.103063 1.000000 -0.065582 -0.004993 -0.005223
X3 0.411351 0.062876 -0.065582 1.000000 -0.038886 -0.178590
X4 0.538308 -0.050210 -0.004993 -0.038886 1.000000 -0.125909
X5 -0.279807 0.061230 -0.005223 -0.178590 -0.125909 1.000000

[Link] 14/30
6/17/25, 6:31 PM gwp1_code

<Figure size 1200x1000 with 0 Axes>

[Link] 15/30
6/17/25, 6:31 PM gwp1_code

Model selection using information criteria (AIC/BIC)

In [29]: print("\n" + "="*50)
print("METHOD 1: MODEL SELECTION USING INFORMATION CRITERIA (AIC/BIC)")
print("="*50)

# List of independent variables

independent_vars = ['X1', 'X2', 'X3', 'X4', 'X5']
best_aic = [Link]
best_bic = [Link]
best_aic_model_name = None
best_bic_model_name = None
best_aic_model = None
best_bic_model = None

all_models_results = []

# Iterate through all possible combinations of independent variables (fro

for i in range(1, len(independent_vars) + 1):
for subset in combinations(independent_vars, i):
# Construct the regression formula
formula = 'Y ~ ' + ' + '.join(subset)

[Link] 16/30
6/17/25, 6:31 PM gwp1_code

# Fit the OLS model

model = [Link](formula, data=data).fit()

# Extract AIC, BIC, and Adjusted R-squared

aic = [Link]
bic = [Link]
adj_r_squared = model.rsquared_adj

model_name = f"Model with {', '.join(subset)}"

all_models_results.append({
'Model': model_name,
'Formula': formula,
'AIC': aic,
'BIC': bic,
'Adjusted R-squared': adj_r_squared,
'Parameters': len(subset) + 1, # +1 for intercept
'Model_Object': model
})

# Update best AIC model

if aic < best_aic:
best_aic = aic
best_aic_model_name = model_name
best_aic_model = model

# Update best BIC model

if bic < best_bic:
best_bic = bic
best_bic_model_name = model_name
best_bic_model = model

# Display Top 5 Models by AIC and BIC

all_models_results_sorted_aic = sorted(all_models_results, key=lambda x:
all_models_results_sorted_bic = sorted(all_models_results, key=lambda x:

print("\nTop 5 Models by AIC:")

for i, res in enumerate(all_models_results_sorted_aic[:5], 1):
print(f"{i}. {res['Model']}: AIC={res['AIC']:.3f}, BIC={res['BIC']:.3

print("\nTop 5 Models by BIC:")

for i, res in enumerate(all_models_results_sorted_bic[:5], 1):
print(f"{i}. {res['Model']}: AIC={res['AIC']:.3f}, BIC={res['BIC']:.3

print(f"\nBest Model by AIC: {best_aic_model_name} (AIC: {best_aic:.3f})"

print(f"Best Model by BIC: {best_bic_model_name} (BIC: {best_bic:.3f})")

# Model summary for the best AIC model

print("\nSummary for Best Model by AIC:")
print(best_aic_model.summary())

# Model summary for the best BIC model (if different from AIC)
if best_aic_model_name != best_bic_model_name:
print("\nSummary for Best Model by BIC:")
print(best_bic_model.summary())

# Create a table of all models with their metrics

model_comparison_df = [Link]([{
'Model': m['Model'],
'Parameters': m['Parameters'],

[Link] 17/30
6/17/25, 6:31 PM gwp1_code

'AIC': m['AIC'],
'BIC': m['BIC'],
'Adj. R-squared': m['Adjusted R-squared']
} for m in all_models_results])

# Sort by AIC
model_comparison_df = model_comparison_df.sort_values('AIC')
print("\nAll Models Comparison (Sorted by AIC):")
print(model_comparison_df.to_string(index=False))

# Visualize AIC and BIC values for all models

[Link](figsize=(12, 6))
[Link](1, 2, 1)
[Link](model_comparison_df['Parameters'], model_comparison_df['AIC']
[Link](model_comparison_df['Parameters'], model_comparison_df['AIC'], '
[Link]('Number of Parameters')
[Link]('AIC')
[Link]('AIC vs Model Complexity')

[Link](1, 2, 2)
[Link](model_comparison_df['Parameters'], model_comparison_df['BIC']
[Link](model_comparison_df['Parameters'], model_comparison_df['BIC'], '
[Link]('Number of Parameters')
[Link]('BIC')
[Link]('BIC vs Model Complexity')

plt.tight_layout()
[Link]()

[Link] 18/30
6/17/25, 6:31 PM gwp1_code

==================================================
METHOD 1: MODEL SELECTION USING INFORMATION CRITERIA (AIC/BIC)
==================================================

Top 5 Models by AIC:

1. Model with X2, X3, X4, X5: AIC=260.617, BIC=273.643, Adj. R-squared=0.6
34
2. Model with X1, X2, X3, X4, X5: AIC=262.593, BIC=278.224, Adj. R-squared
=0.630
3. Model with X2, X3, X4: AIC=264.291, BIC=274.712, Adj. R-squared=0.617
4. Model with X1, X2, X3, X4: AIC=266.191, BIC=279.217, Adj. R-squared=0.6
13
5. Model with X2, X4, X5: AIC=291.890, BIC=302.310, Adj. R-squared=0.495

Top 5 Models by BIC:

1. Model with X2, X3, X4, X5: AIC=260.617, BIC=273.643, Adj. R-squared=0.6
34
2. Model with X2, X3, X4: AIC=264.291, BIC=274.712, Adj. R-squared=0.617
3. Model with X1, X2, X3, X4, X5: AIC=262.593, BIC=278.224, Adj. R-squared
=0.630
4. Model with X1, X2, X3, X4: AIC=266.191, BIC=279.217, Adj. R-squared=0.6
13
5. Model with X2, X4, X5: AIC=291.890, BIC=302.310, Adj. R-squared=0.495

Best Model by AIC: Model with X2, X3, X4, X5 (AIC: 260.617)
Best Model by BIC: Model with X2, X3, X4, X5 (BIC: 273.643)

Summary for Best Model by AIC:

OLS Regression Results
==========================================================================
====
Dep. Variable: Y R-squared:
0.649
Model: OLS Adj. R-squared:
0.634
Method: Least Squares F-statistic: 4
3.87
Date: Tue, 17 Jun 2025 Prob (F-statistic): 8.29
e-21
Time: 18:28:19 Log-Likelihood: -12
5.31
No. Observations: 100 AIC: 2
60.6
Df Residuals: 95 BIC: 2
73.6
Df Model: 4
Covariance Type: nonrobust
==========================================================================
====
coef std err t P>|t| [0.025 0.
975]
--------------------------------------------------------------------------
----
Intercept 1.1893 0.089 13.333 0.000 1.012
1.366
X2 -0.5861 0.091 -6.440 0.000 -0.767 -
0.405
X3 0.5592 0.091 6.124 0.000 0.378
0.740
X4 0.7105 0.082 8.672 0.000 0.548

[Link] 19/30
6/17/25, 6:31 PM gwp1_code

0.873
X5 -0.1966 0.083 -2.355 0.021 -0.362 -
0.031
==========================================================================
====
Omnibus: 4.462 Durbin-Watson:
1.974
Prob(Omnibus): 0.107 Jarque-Bera (JB):
3.920
Skew: 0.473 Prob(JB):
0.141
Kurtosis: 3.215 Cond. No.
1.39
==========================================================================
====

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

All Models Comparison (Sorted by AIC):

Model Parameters AIC BIC Adj. R-sq
uared
Model with X2, X3, X4, X5 5 260.616684 273.642535 0.6
33974
Model with X1, X2, X3, X4, X5 6 262.592528 278.223549 0.6
30170
Model with X2, X3, X4 4 264.291054 274.711735 0.6
16639
Model with X1, X2, X3, X4 5 266.191097 279.216948 0.6
12991
Model with X2, X4, X5 4 291.889597 302.310277 0.4
94796
Model with X1, X2, X4, X5 5 293.801828 306.827678 0.4
89926
Model with X3, X4, X5 4 294.842233 305.262914 0.4
79657
Model with X1, X3, X4, X5 5 296.308890 309.334741 0.4
76977
Model with X3, X4 3 296.442510 304.258021 0.4
66143
Model with X1, X3, X4 4 297.730057 308.150737 0.4
64411
Model with X2, X4 3 299.002018 306.817529 0.4
52303
Model with X1, X2, X4 4 300.988486 311.409167 0.4
46672
Model with X2, X3, X5 4 316.931914 327.352595 0.3
51031
Model with X1, X2, X3, X5 5 318.798302 331.824153 0.3
45076
Model with X4, X5 3 320.380720 328.196231 0.3
21754
Model with X2, X3 3 322.117712 329.933223 0.3
09870
Model with X1, X4, X5 4 322.303799 332.724480 0.3
15216
Model with X1, X2, X3 4 323.812684 334.233364 0.3
04805
Model with X4 2 325.028687 330.239028 0.2

[Link] 20/30
6/17/25, 6:31 PM gwp1_code

82529
Model with X1, X4 3 326.861789 334.677299 0.2
76341
Model with X2, X5 3 331.760939 339.576449 0.2
40005
Model with X1, X2, X5 4 333.760859 344.181540 0.2
32089
Model with X3, X5 3 337.269269 345.084780 0.1
96967
Model with X1, X3, X5 4 338.605294 349.025975 0.1
93972
Model with X2 2 339.907807 345.118147 0.1
67424
Model with X3 2 340.708322 345.918663 0.1
60732
Model with X1, X3 3 341.783826 349.599337 0.1
59883
Model with X1, X2 3 341.867252 349.682762 0.1
59182
Model with X5 2 351.093485 356.303825 0.0
68886
Model with X1, X5 3 352.879853 360.695364 0.0
61295
Model with X1 2 358.868424 364.078765 -0.0
06396

Model selection using LASSO

In [30]: print("\n" + "="*50)
print("METHOD 2: MODEL SELECTION USING LASSO REGRESSION")
print("="*50)

# Prepare data for LASSO

X = data[independent_vars]
y = data['Y']

# Split data into training and testing sets

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2,

# Standardize features
scaler = StandardScaler()

[Link] 21/30
6/17/25, 6:31 PM gwp1_code

X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = [Link](X_test)

# Find the optimal alpha (regularization strength) using cross-validation

lasso_cv = LassoCV(cv=5, random_state=42, max_iter=10000)
lasso_cv.fit(X_train_scaled, y_train)

# Get the best alpha value

best_alpha = lasso_cv.alpha_
print(f"Best Alpha Value: {best_alpha:.6f}")

# Fit LASSO model with the best alpha

lasso_model = Lasso(alpha=best_alpha, random_state=42, max_iter=10000)
lasso_model.fit(X_train_scaled, y_train)

# Get feature coefficients

lasso_coef = [Link]({
'Feature': independent_vars,
'Coefficient': lasso_model.coef_
})

# Sort by absolute coefficient values to identify importance

lasso_coef['Abs_Coefficient'] = [Link](lasso_coef['Coefficient'])
lasso_coef = lasso_coef.sort_values('Abs_Coefficient', ascending=False)

print("\nLASSO Model Coefficients (Sorted by Importance):")

print(lasso_coef)

# Identify selected features (non-zero coefficients)

selected_features = lasso_coef[lasso_coef['Coefficient'] != 0]['Feature']
print(f"\nSelected Features by LASSO: {', '.join(selected_features)}")
print(f"Number of Selected Features: {len(selected_features)} out of {len

# Evaluate the LASSO model on the test set

y_pred_lasso = lasso_model.predict(X_test_scaled)
mse_lasso = mean_squared_error(y_test, y_pred_lasso)
r2_lasso = r2_score(y_test, y_pred_lasso)

print(f"\nLASSO Model Performance:")

print(f"Mean Squared Error: {mse_lasso:.4f}")
print(f"R² Score: {r2_lasso:.4f}")

# Visualize LASSO coefficients

[Link](figsize=(10, 6))
[Link](lasso_coef['Feature'], lasso_coef['Coefficient'])
[Link]('Coefficient Value')
[Link]('Feature')
[Link]('LASSO Regression Coefficients')
[Link](x=0, color='gray', linestyle='--')
[Link](axis='x', linestyle='--', alpha=0.7)
plt.tight_layout()
[Link]()

# Fit LASSO with different alpha values to see how coefficients change
alphas = [Link](-4, 1, 100)
coefs = []

for alpha in alphas:

lasso = Lasso(alpha=alpha, max_iter=10000)
[Link](X_train_scaled, y_train)

[Link] 22/30
6/17/25, 6:31 PM gwp1_code

[Link](lasso.coef_)

# Plot coefficient paths

[Link](figsize=(12, 8))
for i, feature in enumerate(independent_vars):
[Link](alphas, [coef[i] for coef in coefs], label=feature)

[Link]('log')
[Link]('Alpha (Regularization Strength)')
[Link]('Coefficient Value')
[Link]('LASSO Coefficient Paths')
[Link](x=best_alpha, color='k', linestyle='--', label=f'Best Alpha:
[Link](y=0, color='gray', linestyle='--')
[Link]()
[Link](True, alpha=0.3)
plt.tight_layout()
[Link]()

==================================================
METHOD 2: MODEL SELECTION USING LASSO REGRESSION
==================================================
Best Alpha Value: 0.000829

LASSO Model Coefficients (Sorted by Importance):

Feature Coefficient Abs_Coefficient
3 X4 0.789279 0.789279
1 X2 -0.562760 0.562760
2 X3 0.541770 0.541770
4 X5 -0.254143 0.254143
0 X1 -0.082235 0.082235

Selected Features by LASSO: X4, X2, X3, X5, X1

Number of Selected Features: 5 out of 5

LASSO Model Performance:

Mean Squared Error: 0.6833
R² Score: 0.3490

[Link] 23/30
6/17/25, 6:31 PM gwp1_code

Comparison
Key Findings
Information Criteria (AIC/BIC) Approach
Best Model:
Both AIC and BIC selected the same optimal model: Y ~ X2 + X3 + X4 + X5.
Performance Metrics:
Adjusted R²: 0.634 (63.4% of variance explained, highest among candidates).
AIC: 260.617 (lowest, indicating better fit with parsimony).
BIC: 273.643 (lowest, favoring simpler models).
Interpretation:
X1 was excluded as it added negligible explanatory power while increasing
complexity (higher AIC/BIC when included).
The selected model balances goodness-of-fit and parsimony, aligning with CFA
best practices for model selection.

LASSO Regression Approach

Best Model:
Retained all features but shrunk coefficients: Y ~ X4 + X2 + X3 + X5 + X1.
[Link] 24/30
6/17/25, 6:31 PM gwp1_code

Performance Metrics:
R²: 0.349 (weaker than AIC/BIC model).
MSE: 0.6833.
Interpretation:
X1’s coefficient was near zero (-0.082), suggesting minimal impact.
Feature Importance:
X4: Strongest positive effect.
X2: Strongest negative effect.
LASSO’s inclusion of all features (despite shrinkage) hints at potential
multicollinearity or overfitting risks.

Problem 5
5c)
are Why do you
concerned think
about ain Economics
unit root and and
not,Finance
say a we
root
of 1.5?
Simulations to Illustrate the Differences:
Let's simulate both an AR(1) process with a unit root ( ) and an AR(1) process
ϕ = 1

with an explosive root ( ) to visualize their behavior. We'll use standard normal
ϕ = 1.5

white noise ( ).
Wt ∼ N (0, 1)

In [31]: import numpy as np

import [Link] as plt

[Link](42) # For reproducibility

n_steps = 100
white_noise = [Link](0, 1, n_steps)

# Random Walk (phi = 1)

rw_series = [Link](n_steps)
rw_series[0] = 0 # Starting value
for t in range(1, n_steps):
rw_series[t] = rw_series[t-1] + white_noise[t]

# Explosive Process (phi = 1.5)

explosive_series = [Link](n_steps)
explosive_series[0] = 0.1 # Small starting value to show growth
for t in range(1, n_steps):
explosive_series[t] = 1.5 * explosive_series[t-1] + white_noise[t]

[Link](figsize=(12, 6))
[Link](rw_series, label='Random Walk (Unit Root, $\\phi=1$)')
[Link](explosive_series, label='Explosive Process ($\\phi=1.5$)')
[Link]('Simulated Time Series: Unit Root vs. Explosive Process')
[Link]('Time')
[Link]('Value')
[Link]()

[Link] 25/30
6/17/25, 6:31 PM gwp1_code

[Link](True)
[Link]()

Problem 6
In [32]: import numpy as np
import pandas as pd
import [Link] as sm
import [Link] as plt
from scipy import stats

# Set random seed for reproducibility

[Link](42)

# Generate time periods

t = [Link](1, 21)

# Generate X values (explanatory variable)

X = [Link](0, 1, 20)

# Set true parameters

alpha = 2.0
beta1 = 0.5 # Before break
beta2 = 2.0 # After break

# Generate Y with a structural break at t=10

Y = [Link](20)
for i in range(20):
if i < 10: # First period
Y[i] = alpha + beta1 * X[i] + [Link](0, 0.5)
else: # Second period
Y[i] = alpha + beta2 * X[i] + [Link](0, 0.5)

# Create dummy variable for the period after t=10

D = [Link](20)
D[10:] = 1

# Create interaction term

[Link] 26/30
6/17/25, 6:31 PM gwp1_code

X_D = X * D

# Create a DataFrame for better organization

data = [Link]({
't': t,
'Y': Y,
'X': X,
'D': D,
'X_D': X_D
})

print("First few rows of the data:")

print([Link]())
print("\nLast few rows of the data:")
print([Link]())

# Fit the model with the interaction term

X_matrix = sm.add_constant(data[['X', 'X_D']])
model = [Link](Y, X_matrix).fit()

# Print summary statistics

print("\nRegression Results:")
print([Link]())

# Calculate the estimated beta values for each period

beta1_est = [Link][1]
beta2_est = beta1_est + [Link][2]

print(f"\nEstimated parameters:")
print(f"Alpha: {[Link][0]:.4f}")
print(f"Beta1 (t ≤ 10): {beta1_est:.4f}")
print(f"Beta2 (t > 10): {beta2_est:.4f}")
print(f"Change in Beta (Beta2 - Beta1): {[Link][2]:.4f}")

# Test for significance of the structural break

t_stat = [Link][2]
p_value = [Link][2]
print(f"\nTest for structural break:")
print(f"t-statistic for X_D: {t_stat:.4f}")
print(f"p-value: {p_value:.4f}")

if p_value < 0.05:

print("Conclusion: Reject the null hypothesis. There is evidence of a
else:
print("Conclusion: Fail to reject the null hypothesis. No evidence of

# Visualize the results

[Link](figsize=(12, 8))

# Scatter plot of the data points

[Link](X[:10], Y[:10], color='blue', label='Period 1 (t ≤ 10)')
[Link](X[10:], Y[10:], color='red', label='Period 2 (t > 10)')

# Create sequences for plotting the regression lines

X_seq = [Link](min(X) - 0.5, max(X) + 0.5, 100)
Y_seq1 = [Link][0] + [Link][1] * X_seq
Y_seq2 = [Link][0] + [Link][1] * X_seq + [Link][2] * X_

# Plot the regression lines

[Link](X_seq, Y_seq1, 'b-', linewidth=2, label=f'Period 1: Y = {model.p

[Link] 27/30
6/17/25, 6:31 PM gwp1_code

[Link](X_seq, Y_seq2, 'r-', linewidth=2, label=f'Period 2: Y = {model.p

[Link]('Linear Regression with Structural Break at t=10', fontsize=15)

[Link]('X', fontsize=12)
[Link]('Y', fontsize=12)
[Link](True, linestyle='--', alpha=0.7)
[Link](fontsize=12)

# Add a text box with model information

textstr = f'Model: Y(t) = α + β₁X(t) + β₂X_D(t) + ε(t)\n' \
f'α = {[Link][0]:.4f}\n' \
f'β₁ = {beta1_est:.4f}\n' \
f'β₂ = {[Link][2]:.4f}\n' \
f'p-value for β₂ = {p_value:.4f}'

props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)

[Link](0.05, 0.95, textstr, transform=[Link]().transAxes, fontsize=12,
verticalalignment='top', bbox=props)

plt.tight_layout()
[Link]()

[Link] 28/30
6/17/25, 6:31 PM gwp1_code

First few rows of the data:

t Y X D X_D
0 1 2.981181 0.496714 0.0 0.0
1 2 1.817980 -0.138264 0.0 -0.0
2 3 2.357608 0.647689 0.0 0.0
3 4 2.049141 1.523030 0.0 0.0
4 5 1.610732 -0.234153 0.0 -0.0

Last few rows of the data:

t Y X D X_D
15 16 0.265003 -0.562288 1.0 -0.562288
16 17 0.078770 -1.012831 1.0 -1.012831
17 18 1.648660 0.314247 1.0 0.314247
18 19 -0.480141 -0.908024 1.0 -0.908024
19 20 -0.726177 -1.412304 1.0 -1.412304

Regression Results:
OLS Regression Results
==========================================================================
====
Dep. Variable: y R-squared:
0.889
Model: OLS Adj. R-squared:
0.876
Method: Least Squares F-statistic: 6
8.26
Date: Tue, 17 Jun 2025 Prob (F-statistic): 7.52
e-09
Time: 18:28:40 Log-Likelihood: -1
3.083
No. Observations: 20 AIC: 3
2.17
Df Residuals: 17 BIC: 3
5.15
Df Model: 2
Covariance Type: nonrobust
==========================================================================
====
coef std err t P>|t| [0.025 0.
975]
--------------------------------------------------------------------------
----
const 1.8708 0.149 12.577 0.000 1.557
2.185
X 0.3873 0.219 1.771 0.094 -0.074
0.849
X_D 1.5583 0.318 4.894 0.000 0.886
2.230
==========================================================================
====
Omnibus: 1.143 Durbin-Watson:
2.142
Prob(Omnibus): 0.565 Jarque-Bera (JB):
0.678
Skew: 0.445 Prob(JB):
0.713
Kurtosis: 2.854 Cond. No.
4.28
==========================================================================
====

[Link] 29/30
6/17/25, 6:31 PM gwp1_code

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

Estimated parameters:
Alpha: 1.8708
Beta1 (t ≤ 10): 0.3873
Beta2 (t > 10): 1.9456
Change in Beta (Beta2 - Beta1): 1.5583

Test for structural break:

t-statistic for X_D: 4.8938
p-value: 0.0001
Conclusion: Reject the null hypothesis. There is evidence of a structural
break in the beta parameter.

[Link] 30/30