% This program is the base code for testing 
% of Klein's solution algorithm (its complete form). 
%Salemi's Rule

%ORDERING OF THE VARIABLES : yt pt rt yt-1 pt-1 rt-1

clear all;
clc;
close all;
format compact;
warning off;

% Assign values to the model parameters.
nobs=60;
display ('Parameters')



lambda=0.20172729492
a1=1.0957302093
a2=-0.31966056824
b=0.031993 
beta=0.00054 
alfa1=0.605427 
alfa2=0.394903

%Parameters for the loss function

omega_str=zeros(6); %Derived var-cov matrix of the structural shocks
omega_str(1,1)=1.526e-005;
omega_str(1,2)=-1.724e-006;
omega_str(1,3)=6.097e-006;
omega_str(2,1)=-1.724e-006;
omega_str(2,2)=1.384e-005;
omega_str(2,3)=3.418e-006;
omega_str(3,1)=6.097e-006;
omega_str(3,2)=3.418e-006;
omega_str(3,3)=1.778e-005;

delta=.99;

W_tilda=zeros(6,6);
%  W_tilda(1,1)=0.8165; %on y  ln(w_y)=-0.20270824432
%  W_tilda(2,2)=132.7396; %on p ln(w_p)=4.8883895874
%  W_tilda(3,3)=1; %on r
%  W_tilda

LOSS1=zeros(19,18);
LOSS2=zeros(19,18);
LOSS3=zeros(19,18);
LOSS4=zeros(19,18);
LOSS5=zeros(19,18);
LOSS6=zeros(19,18);
LOSS7=zeros(19,18);
comb=0;
Wpindex=1;
Wyindex=1;
for j=0.05:0.05:0.95
  W_tilda(2,2)=j; %on inflation
    for k=0:0.05:(0.95-j)
        W_tilda(1,1)=k; %on output
        W_tilda(3,3)=1-W_tilda(1,1)-W_tilda(2,2); %on interest rate
        W_tilda;
        comb=comb+1

%Find optimal set of thetas and reduced form solution using Klein's algorithm
%by minimizing the loss function.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Rule 1: Full State
theta=[0.419 0.227 0.863 -0.384]; % theta1 (y-1), theta2 (p-1), theta3 (r-1), theta4 (y-2)
display ('FULL STATE: Starting loss, reduced form matrix, eigenvalues,theta, and reduced form and structural covariance matrices');
[loss1, RF1, RF_exp1, roots, M1, L1, N1]=klein1(theta,a1,a2,alfa1,alfa2,b,beta,lambda,omega_str,W_tilda,delta);
theta;
omega1=L1*omega_str(1:3,1:3)*L1';
omega_str;

% Find theta that would minimize the loss function
options=optimset('HessUpdate','bfgs','Display','on','MaxFunEvals'...
    ,1000,'MaxIter',1000,'TolFun',1e-8,'TolX',1e-8,'TolCon',1e-8,'TolPCG',1e-8);
%Other methods of HessUpdate: {'bfgs'} | 'dfp' | 'gillmurray' |'steepdesc'
[theta_opt, fval]=fminunc(@klein1,theta,options,a1,a2,alfa1,alfa2,b,beta,lambda,omega_str,W_tilda,delta);

display ('Optimal loss, reduced form, eigenvalues, and thetas');
[loss1, RF1, RF_exp1, roots, M1, L1, N1]=klein1(theta_opt,a1,a2,alfa1,alfa2,b,beta,lambda,omega_str,W_tilda,delta);
theta_opt;
omega1=L1*omega_str(1:3,1:3)*L1';

LOSS1(Wyindex,Wpindex)=loss1;

data1(comb,1)=W_tilda(2,2); 
data1(comb,2)=W_tilda(1,1);
data1(comb,3)=W_tilda(3,3);
data1(comb,4)=loss1;
data1(comb,5)=M1(2,2);
data1(comb,6)=M1(1,1);
data1(comb,7)=M1(3,3);
data1(comb,8)=theta_opt(2); 
data1(comb,9)=theta_opt(1);
data1(comb,10)=theta_opt(3);
data1(comb,11)=theta_opt(4); 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Rule 2: Taylor (backward-looking)
theta=[0.4 0.5]; % theta1 (y-1), theta2 (p-1)
display ('TAYLOR BACKWARD-LOOKING: Starting loss, reduced form matrix, eigenvalues,theta, and reduced form and structural covariance matrices');
[loss2, RF2, RF_exp2, roots, M2, L2, N2]=klein2(theta,a1,a2,alfa1,alfa2,b,beta,lambda,omega_str,W_tilda,delta);
theta;
omega2=L2*omega_str(1:3,1:3)*L2';
omega_str;

% Find theta that would minimize the loss function
options=optimset('HessUpdate','bfgs','Display','on','MaxFunEvals'...
    ,1000,'MaxIter',1000,'TolFun',1e-8,'TolX',1e-8,'TolCon',1e-8,'TolPCG',1e-8);
%Other methods of HessUpdate: {'bfgs'} | 'dfp' | 'gillmurray' |'steepdesc'
[theta_opt, fval]=fminunc(@klein2,theta,options,a1,a2,alfa1,alfa2,b,beta,lambda,omega_str,W_tilda,delta);

display ('Optimal loss, reduced form, eigenvalues, and thetas');
[loss2, RF2, RF_exp2, roots, M2, L2, N2]=klein2(theta_opt,a1,a2,alfa1,alfa2,b,beta,lambda,omega_str,W_tilda,delta);
theta_opt;
omega2=L2*omega_str(1:3,1:3)*L2';

LOSS2(Wyindex,Wpindex)=loss2;

data2(comb,1)=W_tilda(2,2);
data2(comb,2)=W_tilda(1,1);
data2(comb,3)=W_tilda(3,3);
data2(comb,4)=loss2;
data2(comb,5)=M2(2,2);
data2(comb,6)=M2(1,1);
data2(comb,7)=M2(3,3);
data2(comb,8)=theta_opt(2);
data2(comb,9)=theta_opt(1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Rule 3: Optimal Taylor (current y and p)
theta=[0.4 0.5] ;% theta1 (y), theta2 (p)
display ('OPTIMAL TAYLOR (CURRENT Y AND P): Starting loss, reduced form matrix, eigenvalues,theta, and reduced form and structural covariance matrices');
[loss3, RF3, RF_exp3, roots, M3, L3, N3]=klein3(theta,a1,a2,alfa1,alfa2,b,beta,lambda,omega_str,W_tilda,delta);
theta;
omega3=L3*omega_str(1:3,1:3)*L3';
omega_str;

% Find theta that would minimize the loss function
options=optimset('HessUpdate','bfgs','Display','on','MaxFunEvals'...
    ,1000,'MaxIter',1000,'TolFun',1e-8,'TolX',1e-8,'TolCon',1e-8,'TolPCG',1e-8);
%Other methods of HessUpdate: {'bfgs'} | 'dfp' | 'gillmurray' |'steepdesc'
[theta_opt, fval]=fminunc(@klein3,theta,options,a1,a2,alfa1,alfa2,b,beta,lambda,omega_str,W_tilda,delta);

display ('Optimal loss, reduced form, eigenvalues, and thetas');
[loss3, RF3, RF_exp3, roots, M3, L3, N3]=klein3(theta_opt,a1,a2,alfa1,alfa2,b,beta,lambda,omega_str,W_tilda,delta);
theta_opt;
omega3=L3*omega_str(1:3,1:3)*L3';

LOSS3(Wyindex,Wpindex)=loss3;

data3(comb,1)=W_tilda(2,2);
data3(comb,2)=W_tilda(1,1);
data3(comb,3)=W_tilda(3,3);
data3(comb,4)=loss3;
data3(comb,5)=M3(2,2);
data3(comb,6)=M3(1,1);
data3(comb,7)=M3(3,3);
data3(comb,8)=theta_opt(2);
data3(comb,9)=theta_opt(1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Rule 4: Optimal Taylor with Interest Rate Smooting (current y and p)
theta=[0.2 0.4 0.3]; % theta1 (y), theta2 (p), theta3 (r-1)
display ('OPTIMAL TAYLOR WITH INTEREST RATE SMOOTHING (CURRENT Y AND P): Starting loss, reduced form matrix, eigenvalues,theta, and reduced form and structural covariance matrices');
[loss4, RF4, RF_exp4, roots, M4, L4, N4]=klein4(theta,a1,a2,alfa1,alfa2,b,beta,lambda,omega_str,W_tilda,delta);
theta;
omega4=L4*omega_str(1:3,1:3)*L4';
omega_str;

% Find theta that would minimize the loss function
options=optimset('HessUpdate','bfgs','Display','on','MaxFunEvals'...
    ,1000,'MaxIter',1000,'TolFun',1e-8,'TolX',1e-8,'TolCon',1e-8,'TolPCG',1e-8);
%Other methods of HessUpdate: {'bfgs'} | 'dfp' | 'gillmurray' |'steepdesc'
[theta_opt, fval]=fminunc(@klein4,theta,options,a1,a2,alfa1,alfa2,b,beta,lambda,omega_str,W_tilda,delta);

display ('Optimal loss, reduced form, eigenvalues, and thetas');
[loss4, RF4, RF_exp4, roots, M4, L4, N4]=klein4(theta_opt,a1,a2,alfa1,alfa2,b,beta,lambda,omega_str,W_tilda,delta);
theta_opt;
omega4=L4*omega_str(1:3,1:3)*L4';

LOSS4(Wyindex,Wpindex)=loss4;

data4(comb,1)=W_tilda(2,2);
data4(comb,2)=W_tilda(1,1);
data4(comb,3)=W_tilda(3,3);
data4(comb,4)=loss4;
data4(comb,5)=M4(2,2);
data4(comb,6)=M4(1,1);
data4(comb,7)=M4(3,3);
data4(comb,8)=theta_opt(2);
data4(comb,9)=theta_opt(1);
data4(comb,10)=theta_opt(3);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Rule 5: Taylor with the Original Coefficients (current y and p)
theta=[0.5 1.5]; % theta1 (y), theta2 (p)
display ('TAYLOR WITH THE ORIGINAL COEFFICIENTS (CURRENT Y AND P): Starting loss, reduced form matrix, eigenvalues,theta, and reduced form and structural covariance matrices');
[loss5, RF5, RF_exp5, roots, M5, L5, N5]=klein5(theta,a1,a2,alfa1,alfa2,b,beta,lambda,omega_str,W_tilda,delta);
omega5=L5*omega_str(1:3,1:3)*L5';
omega_str;

LOSS5(Wyindex,Wpindex)=loss5;

data5(comb,1)=W_tilda(2,2);
data5(comb,2)=W_tilda(1,1);
data5(comb,3)=W_tilda(3,3);
data5(comb,4)=loss5;
data5(comb,5)=M5(2,2);
data5(comb,6)=M5(1,1);
data5(comb,7)=M5(3,3);
data5(comb,8)=theta(2);
data5(comb,9)=theta(1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Rule 6: Goodhart (current p)
theta=[0.4]; % theta1 (p)
display ('GOODHART (CURRENT P): Starting loss, reduced form matrix, eigenvalues,theta, and reduced form and structural covariance matrices');
[loss6, RF6, RF_exp6, roots, M6, L6, N6]=klein6(theta,a1,a2,alfa1,alfa2,b,beta,lambda,omega_str,W_tilda,delta);
theta;
omega6=L6*omega_str(1:3,1:3)*L6';
omega_str;

% Find theta that would minimize the loss function
options=optimset('HessUpdate','bfgs','Display','on','MaxFunEvals'...
    ,1000,'MaxIter',1000,'TolFun',1e-8,'TolX',1e-8,'TolCon',1e-8,'TolPCG',1e-8);
%Other methods of HessUpdate: {'bfgs'} | 'dfp' | 'gillmurray' |'steepdesc'
[theta_opt, fval]=fminunc(@klein6,theta,options,a1,a2,alfa1,alfa2,b,beta,lambda,omega_str,W_tilda,delta);

display ('Optimal loss, reduced form, eigenvalues, and thetas');
[loss6, RF6, RF_exp6, roots, M6, L6, N6]=klein6(theta_opt,a1,a2,alfa1,alfa2,b,beta,lambda,omega_str,W_tilda,delta);
theta_opt;
omega6=L6*omega_str(1:3,1:3)*L6';

LOSS6(Wyindex,Wpindex)=loss6;

data6(comb,1)=W_tilda(2,2);
data6(comb,2)=W_tilda(1,1);
data6(comb,3)=W_tilda(3,3);
data6(comb,4)=loss6;
data6(comb,5)=M6(2,2);
data6(comb,6)=M6(1,1);
data6(comb,7)=M6(3,3);
data6(comb,8)=theta_opt(1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Rule 7: Goodhart (forward-looking p)
theta=[0.9] ;% theta1 (p+1)
display ('GOODHART (FORWARD-LOOKING P): Starting loss, reduced form matrix, eigenvalues,theta, and reduced form and structural covariance matrices');
[loss7, RF7, RF_exp7, roots, M7, L7, N7]=klein7(theta,a1,a2,alfa1,alfa2,b,beta,lambda,omega_str,W_tilda,delta);
theta;
omega7=L7*omega_str(1:3,1:3)*L7';
omega_str;

% Find theta that would minimize the loss function
options=optimset('HessUpdate','bfgs','Display','on','MaxFunEvals'...
    ,1000,'MaxIter',1000,'TolFun',1e-8,'TolX',1e-8,'TolCon',1e-8,'TolPCG',1e-8);
%Other methods of HessUpdate: {'bfgs'} | 'dfp' | 'gillmurray' |'steepdesc'
[theta_opt, fval]=fminunc(@klein7,theta,options,a1,a2,alfa1,alfa2,b,beta,lambda,omega_str,W_tilda,delta);

display ('Optimal loss, reduced form, eigenvalues, and thetas');
[loss7, RF7, RF_exp7, roots, M7, L7, N7]=klein7(theta_opt,a1,a2,alfa1,alfa2,b,beta,lambda,omega_str,W_tilda,delta);
theta_opt;
omega7=L7*omega_str(1:3,1:3)*L7';

LOSS7(Wyindex,Wpindex)=loss7;

data7(comb,1)=W_tilda(2,2);
data7(comb,2)=W_tilda(1,1);
data7(comb,3)=W_tilda(3,3);
data7(comb,4)=loss7;
data7(comb,5)=M7(2,2);
data7(comb,6)=M7(1,1);
data7(comb,7)=M7(3,3);
data7(comb,8)=theta_opt(1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Wyindex=Wyindex+1;
   end; % loop on output
Wyindex=1;
%saving data to file
save rule1 data1 LOSS1 -ASCII;
save rule2 data2 LOSS2 -ASCII;
save rule3 data3 LOSS3 -ASCII;
save rule4 data4 LOSS4 -ASCII;
save rule5 data5 LOSS5 -ASCII;
save rule6 data6 LOSS6 -ASCII;
save rule7 data7 LOSS7 -ASCII;
   
Wpindex=Wpindex+1;   
end; %loop on inflation

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Building impulse response functions

X=zeros(4,nobs); % Reduced form sollution for predetermined variables
EXP=zeros(2,nobs); % Reduced form solution for NOT predetermined variables
X_lag=zeros(4,nobs);

str_shock=[0.01; 0.0; 0.0; 0.0]; % Structural shocks (on y, p, r, lagged y)

for i=2:nobs
    %Computing impulse response functions for the predetermined variables
    X(:,i)=real(RF1(1:4,1:4))*X_lag(:,i-1)+real(L1)*str_shock(1:3); 
    %Computing impulse response functions for the NOT predetermined variables
    EXP(:,i)=real(RF_exp1)*X_lag(:,i-1)+real(N1)*str_shock(1:3); %expected y, expected p

    X_lag(1,i)=X(1,i); %lagged y
    X_lag(2,i)=X(2,i); %lagged p
    X_lag(3,i)=X(3,i);%lagged r
    X_lag(4,i)=X(4,i-1); %twice lagged y
    
    str_shock=zeros(4,1);
end;


%Impulse responce functions for predetermined variables:
figure
subplot (3,1,1)
plot ([X(1,:)' zeros(nobs,1)]) 
title ('Output')
subplot (3,1,2)
plot ([X(2,:)' zeros(nobs,1)])
title ('Inflation')
subplot (3,1,3)
plot ([X(3,:)' zeros(nobs,1)])
title ('Interest rate')

%Impulse responce functions for NOT predetermined variables:
figure
subplot (2,1,1)
plot ([EXP(1,:)' zeros(nobs,1)]) 
title ('Expected Output')
subplot (2,1,2)
plot ([EXP(2,:)' zeros(nobs,1)])
title ('Expected Inflation')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Choosing the best

all=zeros(19,18,7);
all(:,:,1)=LOSS1;
all(:,:,2)=LOSS2;
all(:,:,3)=LOSS3; 
all(:,:,4)=LOSS4;
all(:,:,5)=LOSS5;
all(:,:,6)=LOSS6;
all(:,:,7)=LOSS7; 

[best, index]=min(all,[],3) %best - the minimum loss, index - the rule under which it can be obtained
save best best index -ascii;