%% Parameters calibration (Gali ch. 3 pag. 52)

alpha = 1/3; 
sigma = 1; % log utility
phi = 1; % unitary Frisch elasticity
epsilon = 6;
theta = 2/3; % average price duration of three quarters
beta = 0.99; 
rho = 0.5; % shock's persistence
theta_cap = (1 - alpha) / (1 - alpha + alpha*epsilon);
lambda = theta_cap * ((1 - theta) * ( 1 - beta*theta) / theta);
k = lambda * (sigma + (phi + alpha) / (1 - alpha));
alpha_x = k / epsilon;
rho_mp = 0.5; % mp shock persistence
psi_cp = (1/(alpha_x*(1-beta*rho)+k^2));
psi_mp = (1/(alpha_x*(1-beta*rho_mp)+k^2));

% load stored random numbers 
load('cost_shock.mat')
load('mp_shock.mat')

cost_shock_ar1(1) = cost_shock(1);
for i = 1:999
    cost_shock_ar1(i+1) = rho*cost_shock_ar1(i) + cost_shock(i);
end

mp_shock_ar1(1) = mp_shock(1);

for i=1:999
mp_shock_ar1(i+1)=rho_mp*mp_shock_ar1(i) + mp_shock(i);
end

%
%% Optimal Monetary Policy Under Discretion
%% 1.  Cost-push shock only

pi = (alpha_x / (alpha_x * (1-rho*beta) + k ^ 2)) * cost_shock_ar1;
x = - (k / alpha_x) * pi;

% Chart
chart1 = figure;
plot(1:10);
scatter(x, pi, 40, 'ro');
xlabel ('Output gap (%)')
ylabel ('Inflation (percentage point deviation from target)')
z = lsline
set(z(1),'color','r', 'LineWidth', 2.5)
axis([-10 10 -2.5 2.5])
lsline
set(chart1, 'Units', 'Inches');
pos = get(chart1,'Position');
set(chart1, 'PaperPositionMode', 'Auto', 'PaperUnits', 'Inches','Papersize', [pos(3),pos(4)])
print(chart1, 'Chart1', '-dpdf', '-r0')

% Regressions
beta_hat = fitlm(x, pi, 'Intercept', false)

%
%% 2. Monetary policy and cost-push shocks NOT controlling for expectations (NEED TO REVIEW EXP PI BOTH

% Case 1: equal relative variance
pi_both = alpha_x*psi_cp*cost_shock_ar1 + (k^2)*psi_mp*mp_shock_ar1;
x_both = -k*psi_cp*cost_shock_ar1 + (k/alpha_x)*(1-(k^2)*psi_mp)*mp_shock_ar1;

exp_pi_both = ((k^2)/alpha_x)*((alpha_x + k^2)*psi_mp -1)*mp_shock_ar1 + ((alpha_x+k^2)*psi_cp - 1)*cost_shock_ar1;
X1 = [x_both, exp_pi_both];

modelX = fitlm(X1,pi_both, 'Intercept', false)

% Case 2: higher relative variance of the cost-push shock 
mp_shock_small = (1/100) .* mp_shock_ar1;
pi_both_highvarcp = alpha_x*psi_cp*cost_shock_ar1 + (k^2)*psi_mp*mp_shock_small;
x_both_highvarcp = -k*psi_cp*cost_shock_ar1 + (k/alpha_x)*(1-(k^2)*psi_mp)*mp_shock_small;



% Case 3: higher relative variance of mp shock
cost_shock_small = (1/100).*cost_shock_ar1;
pi_both_highvarmp = alpha_x*psi_cp*cost_shock_small + (k^2)*psi_mp*mp_shock_ar1;
x_both_highvarmp = -k*psi_cp*cost_shock_small + (k/alpha_x)*(1-(k^2)*psi_mp)*mp_shock_ar1;

exp_pi_both_highvarmp = ((k^2)/alpha_x)*((alpha_x + k^2)*psi_mp -1)*mp_shock_ar1 + ((alpha_x+k^2)*psi_cp - 1)*cost_shock_small;
X3 = [x_both_highvarmp, exp_pi_both_highvarmp];

modelX3 = fitlm(X3,pi_both_highvarmp, 'Intercept', false)


% Chart
chart2 = figure
plot(1:10);
scatter(x_both_highvarcp,pi_both_highvarcp, 40,'ro')
lsline
hold on
scatter(x_both_highvarmp,pi_both_highvarmp, 40, 'bo')
lsline
hold on
scatter(x_both, pi_both, 40, 'go')
lsline
axis([-10 10 -2.5 2.5])    
h = lsline
set(h(1),'color','g', 'LineWidth', 2.5)
set(h(2),'color','b', 'LineWidth', 2.5)
set(h(3),'color','r', 'LineWidth', 2.5)
ylabel('Inflation (percentage point deviation from target)')
xlabel('Output gap (%)')
set(chart2, 'Units', 'Inches');
pos2 = get(chart2, 'Position');
set(chart2, 'PaperPositionMode', 'Auto', 'PaperUnits', 'Inches', 'PaperSize' , [pos2(3),pos(4)])
print(chart2, 'Chart2', '-dpdf', '-r0')

% Regressions
beta_both = fitlm(x_both, pi_both, 'Intercept', false)
beta_both_highvarcp = fitlm(x_both_highvarcp, pi_both_highvarcp, 'Intercept', false)
beta_both_highvarmp = fitlm(x_both_highvarmp, pi_both_highvarmp, 'Intercept', false)

%
%% 3. Monetary policy and cost-push shocks CONTROLLING for expectations

% Case 1: equal relative variance
pi_minus_exp_both = (k^2)*(1-beta*rho_mp)*psi_mp*mp_shock_ar1 - (k^2)*psi_cp*(1-(1/psi_cp))*cost_shock_ar1;
x_minus_exp_both = k*(1-beta*rho_mp)*psi_mp*mp_shock_ar1 - (k*psi_cp*(1-(1/psi_cp))-(1/k))*cost_shock_ar1;


% Case 2: higher relative variance of the cost-push shock
pi_minus_exp_both_highvarcp = (k^2)*(1-beta*rho_mp)*psi_mp*mp_shock_small - (k^2)*psi_cp*(1-(1/psi_cp))*cost_shock_ar1;
x_minus_exp_both_highvarcp = k*(1-beta*rho_mp)*psi_mp*mp_shock_small - (k*psi_cp*(1-(1/psi_cp))-(1/k))*cost_shock_ar1;

% Case 3: higher relative variance of the mp shock
pi_minus_exp_both_highvarmp = (k^2)*(1-beta*rho_mp)*psi_mp*mp_shock_ar1 - (k^2)*psi_cp*(1-(1/psi_cp))*cost_shock_small;
x_minus_exp_both_highvarmp = k*(1-beta*rho_mp)*psi_mp*mp_shock_ar1 - (k*psi_cp*(1-(1/psi_cp))-(1/k))*cost_shock_small;

% Regressions
beta_minus_exp_both = fitlm(x_minus_exp_both, pi_minus_exp_both, 'Intercept', false)
beta_minus_exp_both_highvarcp = fitlm(x_minus_exp_both_highvarcp, pi_minus_exp_both_highvarcp, 'Intercept', false)
beta_minus_exp_both_highvarmp = fitlm(x_minus_exp_both_highvarmp, pi_minus_exp_both_highvarmp, 'Intercept', false)

%
%% Optimal Monetary Policy Under Committment
%
%% Extra parameters

a = alpha_x / (alpha_x*(1+beta) + k^2);
delta = ((1 - sqrt(1 - 4*beta*a^2)) / (2*a*beta));

%
%% Equilibrium price level and output gap - pag. 104 eq (15) & (16)

price_level(1) = (delta / (1 - delta*beta*rho) * cost_shock_ar1(1))';
for i=1:999
    price_level(i+1) = delta * price_level(i) + (delta / (1 - delta*beta*rho) * cost_shock_ar1(i+1))';
end

x_commitment(1) = - (k*delta / (alpha_x*(1-delta*beta*rho))) * cost_shock_ar1(1);
for i=1:999
    x_commitment(i+1) = delta*x(i) - (k*delta / (alpha_x*(1-delta*beta*rho))) * cost_shock_ar1(i+1);
end



%
%% Regression

p_1 = [0;price_level(1:end-1)']; % shifting price level down by one element
pi_comm = [price_level' - p_1]; % calculating price level difference
pi_comm2 = pi_comm(2:end); % removing the first value from the matrix
x_comm = x_commitment(2:end)'; % removing first value from the output gap matrix

beta_hat_comm = (x_comm' * x_comm)^(-1) * x_comm' * pi_comm2

chart3 = figure
scatter(x_comm,pi_comm2, 'ro')
xlabel('Output gap (%)')
ylabel('Inflation rate (percentage point deviations from target)')
lsline
f = lsline
set(f(1),'color','r', 'LineWidth', 2.5)
axis([-10 10 -2.5 2.5])
lsline
set(chart3, 'Units', 'Inches');
pos = get(chart3,'Position');
set(chart3, 'PaperPositionMode', 'Auto', 'PaperUnits', 'Inches','Papersize', [pos(3),pos(4)])
print(chart3, 'Chart3', '-dpdf', '-r0')


figure
scatter(x, price_level,'y.')
xlabel('Output gap (% deviations from SS output gap)')
ylabel('Price level (deviations from SS price level)')