clear all;
close all;
clc;

M = importdata('data2.xlsx',',',1);
data = M.data;
n = size(data,1);

%% Rename Variables 
time = data(:,1);
%C = b1 + b2*Y
x = data(:,4); %income
y = data(:,3); %consumption

constant=ones(n,1); %vector of 1s

xmat=[constant,x]; %the x matrix
[n,k] = size(xmat); %size


scatter(x, y)
ylabel({'Consumption'});
xlabel({'Income'});



%% Standard Error of MPC

b = inv(xmat' * xmat) * (xmat' * y);
disp('The autonomous consumption is = '); disp(b(1)) %Display intercept
disp('The MPC is = '); disp(b(2))%Display slope


e = y - xmat*b; %error
s2 = (e'*e)/(n-k); %sigma square
sigma = sqrt((e'*e)/(n-k));
varb = s2*inv(xmat'*xmat); %variance produces a variance-covariance matrix 
stder=sqrt(varb); %Taking square root to get the standard error
SE=diag(stder); 
disp('The sigma is = '); disp(sigma);
%disp('The Standard errors is = '); disp(SE);

%% The govt multiplier
gm = 1./(1-b(2));
disp('The Government Multiplier is = '); disp(gm)

% We cant use the same formula to find the var(gm) 
% Hence we use the first-order Taylor series expansion


% % Loop to Simulate Data 
m= 1000;
b_hat = zeros(m,2);
gamma_hat = zeros(m,1);
for i = (1:1:m);
    y = xmat*b + sigma.*randn(n,1);
    b_hat(i,:) = (inv(xmat'*xmat)*(xmat'*y))';
    gamma_hat(i,:) = 1./(1-b_hat(i,2));
end;

se_mc_gamma_hat = sqrt(var(gamma_hat));

% % Graph Sampling Distributions
 hist(gamma_hat, 100)

% %delta: se
v_gamma_hat = ((1-b(2)).^(-4))*varb(2,2);
se_gamma_hat = sqrt(((1-b(2)).^(-4))*varb(2,2));
% 
disp('The variance of Gamma from Delta Method is = '); disp(v_gamma_hat);
disp('The Standard error of Gamma from the Delta Method is = '); disp(se_gamma_hat);
%