%% Week 11: Empirical Macro - OLS and VARs
% This script covers:
% 1. Data loading and transformation
% 2. OLS estimation (Phillips Curve)
% 3. VAR estimation (Bivariate)
% 4. Impulse Response Functions (IRFs)

clear; clc; close all;

if ~exist('Figures','dir'), mkdir Figures; end


%% 1. Load Data
% Ensure you have run week11_get_data.m first to generate the CSV.
if ~isfile('Data/US_macro_data.csv')
    error('Data/US_macro_data.csv not found. Run week11_get_data.m first.');
end

T = readtable('Data/US_macro_data.csv');
T.Observation_date = datetime(T.Observation_date);

% Extract series
gdp_level = T.GDPC1;
deflator  = T.GDPDEF;
ffr       = T.FEDFUNDS;
dates     = T.Observation_date;

%% 2. Transformations
% We need stationary series for VARs.
% GDP -> Log Difference (Growth Rate)
% Deflator -> Log Difference (Inflation)
% Interest Rate -> Levels (already percent)

% Annualized Quarterly Growth Rate: 400 * (ln(Y_t) - ln(Y_{t-1}))
gdp_growth = [NaN; 400 * diff(log(gdp_level))];

% Annualized Inflation Rate: 400 * (ln(P_t) - ln(P_{t-1}))
inflation  = [NaN; 400 * diff(log(deflator))];

% Align data (remove the first NaN observation)
% Data Matrix: [GDP_Growth, Inflation, FFR]
data_raw = [gdp_growth, inflation, ffr];

dates_clean = dates(2:end);
data_clean  = data_raw(2:end, :);

% Plot transformed data
gfc=figure('Name', 'Transformed Data');
subplot(3,1,1); plot(dates_clean, data_clean(:,1)); title('Real GDP Growth (Ann %)'); grid on;
subplot(3,1,2); plot(dates_clean, data_clean(:,2)); title('Inflation (GDP Deflator, Ann %)'); grid on;
subplot(3,1,3); plot(dates_clean, data_clean(:,3)); title('Federal Funds Rate (%)'); grid on;

exportgraphics(gfc,'Figures/trasformed_data.png')

%% 3. OLS: The Phillips Curve
% Simple specification: pi_t = alpha + beta * y_t + e_t
% (Note: In reality, we'd use output gap, but growth is a proxy for this exercise)

% WARNING: This specification likely suffers from ENDOGENEITY (Simultaneity Bias).
% Inflation and GDP are jointly determined. OLS estimates here are illustrative
% of "naive" correlations, motivating the need for VARs and structural identification.

Y_ols = data_clean(:, 2); % Inflation
X_ols = [ones(length(Y_ols),1), data_clean(:, 1)]; % Constant + GDP Growth

% Manual OLS Estimator: beta_hat = (X'X)^-1 X'Y
beta_hat = (X_ols' * X_ols) \ (X_ols' * Y_ols);

fprintf('\n--- OLS Results (Manual) ---\n');
fprintf('Constant: %.4f\n', beta_hat(1));
fprintf('Slope (GDP Growth): %.4f\n', beta_hat(2));

% Compare with MATLAB built-in fitlm
mdl = fitlm(data_clean(:,1), data_clean(:,2));
disp(mdl);

% Plot scatter and regression line
gfc1=figure('Name', 'Phillips Curve OLS');
scatter(data_clean(:,1), data_clean(:,2), 'filled'); hold on;
plot(data_clean(:,1), X_ols * beta_hat, 'r-', 'LineWidth', 2);
xlabel('GDP Growth'); ylabel('Inflation');
title('Simple Phillips Curve Relationship');
grid on;
exportgraphics(gfc1,'Figures/PC.png')

%% 4. Bivariate VAR(1) Estimation
% System: Y_t = [GDP_Growth_t, Inflation_t]'
% Model: Y_t = A0 + A1 * Y_{t-1} + e_t

% Create Lagged Matrices
% Y_var (Dependent): t = 2 to T
% X_var (Regressors): Constant + Y_{t-1} (t = 1 to T-1)

Y_vec = data_clean(:, 1:2); % Use only GDP and Inflation
T_sample = size(Y_vec, 1);

Y_lhs = Y_vec(2:end, :);      % t = 2...T
X_rhs = [ones(T_sample-1, 1), Y_vec(1:end-1, :)]; % Constant, t = 1...T-1

% Estimate coefficients equation by equation
% Equation 1: GDP Growth
B_gdp = (X_rhs' * X_rhs) \ (X_rhs' * Y_lhs(:,1));

% Equation 2: Inflation
B_inf = (X_rhs' * X_rhs) \ (X_rhs' * Y_lhs(:,2));

% Display Coefficients
fprintf('\n--- VAR(1) Results ---\n');
fprintf('GDP Eq:  c=%.2f, b_gdp=%.2f, b_inf=%.2f\n', B_gdp(1), B_gdp(2), B_gdp(3));
fprintf('Inf Eq:  c=%.2f, b_gdp=%.2f, b_inf=%.2f\n', B_inf(1), B_inf(2), B_inf(3));

% Construct Matrix A1 for dynamics
% Y_t = A0 + A1 * Y_{t-1}
A1 = [B_gdp(2:3)'; B_inf(2:3)']; % Transpose to get row per equation
disp('Transition Matrix A1:');
disp(A1);

%% 5. Residuals and Structural Identification
% Calculate Residuals
res_gdp = Y_lhs(:,1) - X_rhs * B_gdp;
res_inf = Y_lhs(:,2) - X_rhs * B_inf;
residuals = [res_gdp, res_inf];

% Covariance Matrix (Reduced Form Sigma)
Sigma = cov(residuals);
disp('Residual Covariance Matrix (Sigma):');
disp(Sigma);

% Cholesky Decomposition: Sigma = P * P'
% P is Lower Triangular.
% Ordering: GDP first, Inflation second.
% Interpretation:
% - GDP shock affects GDP immediately, but Inflation only with lag (via VAR).
% - Inflation shock affects Inflation immediately, AND GDP immediately (Wait, check lower triangular def)
% Lower Triangular P:
% [ p11  0  ]  [ u1 ]  -> y1 shock affects only y1 at t=0
% [ p21 p22 ]  [ u2 ]  -> y2 shock affects y2, AND y1 shock affects y2 at t=0
%
% Standard Cholesky implies:
% Variable 1 (GDP) does NOT respond to Variable 2 (Inf) at t=0.
% Variable 2 (Inf) DOES respond to Variable 1 (GDP) at t=0.
P = chol(Sigma, 'lower');
disp('Cholesky Factor P (Lower):');
disp(P);

%% 6. Impulse Response Functions (IRFs)
H = 16; % Horizon (quarters)
n_vars = 2;

% Storage for IRFs
% We will compute responses to:
% Shock 1: GDP Shock (1 std dev) -> Vector [P(1,1); P(2,1)]
% Shock 2: Inflation Shock (1 std dev) -> Vector [0; P(2,2)] (Strictly, column 2 of P)

shock_1 = P(:, 1); % Structural GDP shock
shock_2 = P(:, 2); % Structural Inflation shock

% Initialize IRF paths
irf_1 = zeros(n_vars, H);
irf_2 = zeros(n_vars, H);

% Impact (t=1)
irf_1(:, 1) = shock_1;
irf_2(:, 1) = shock_2;

% Iterate forward: x_{t+1} = A1 * x_t
for h = 2:H
    irf_1(:, h) = A1 * irf_1(:, h-1);
    irf_2(:, h) = A1 * irf_2(:, h-1);
end

%% 7. Plot IRFs
t_axis = 0:H-1;

gfc2=figure('Name', 'VAR Impulse Responses');

% Response to GDP Shock
subplot(2,2,1); plot(t_axis, irf_1(1,:), 'b-o', 'LineWidth', 1.5); 
yline(0, 'k--'); title('Response of GDP to GDP Shock'); grid on; xlim([0 H-1]);
subplot(2,2,2); plot(t_axis, irf_1(2,:), 'r-o', 'LineWidth', 1.5); 
yline(0, 'k--'); title('Response of Inflation to GDP Shock'); grid on; xlim([0 H-1]);

% Response to Inflation Shock
subplot(2,2,3); plot(t_axis, irf_2(1,:), 'b-o', 'LineWidth', 1.5); 
yline(0, 'k--'); title('Response of GDP to Inf Shock'); grid on; xlim([0 H-1]);
subplot(2,2,4); plot(t_axis, irf_2(2,:), 'r-o', 'LineWidth', 1.5); 
yline(0, 'k--'); title('Response of Inflation to Inf Shock'); grid on; xlim([0 H-1]);

exportgraphics(gfc2,'Figures/VAR2.png')