%Program to solve linear model with higher order expectation of current and future endogenous variables (c) K Nimark

%change parameters in lines 7 - 15 to experiment with different calibrations
 clear all;
% Choose kbar and tolerance
kbar=15; %Choose number of orders of expectations to be included
tol=0.00000001;Step=.9;
%Define and assign parameter values%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
theta=.9; beta=0.995;rho=.9;gamma=2;varphi=2; fipi=1.5; fiy=.5; %fiy=0 does not substantially change results
      %variances
      sigma2e=1; %Too large value leads to matrix inversion problems
      sigma2v=1;  
      sigma2y=1;  
      sigma2pi=0.01;%Too small value leads to matrix inversion problems
      SigmaEE = diag([sigma2e,sigma2pi,sigma2y]);
      SigmaUU = diag([sigma2v,sigma2pi,sigma2y]);
      
for i=1:kbar;
    a(1,i)=((1-theta)*(1-theta*beta))*((1-theta)^(i-1));
    b(1,i)=(beta*theta)*(1-theta)^(i-1);    
end;

%Assign starting values to M,N,c,d, D and SigmaVV%%%%%%%%%%%%%%%%%%%%%%%%%%
M=zeros(2*kbar,2*kbar);M(1,1)=rho;
Mst=M;
N=zeros(4*kbar,3);N(1:2,1:2)=eye(2);

c=zeros(1,kbar*2);c(1,1)=a(1,1);c(1,2)=a(1,1);
cst=c;

dst=zeros(1,kbar*2);
ddiff=1;
while ddiff > tol;
    d=dst*M-(fipi/gamma)*c+(1/gamma)*c*M;
    ddiff=max(abs(d-dst));
    dst= Step*d+ (1-Step)*dst;
end;
d=dst;

%Construct G matrix (in the text G is the part of the matrix D in (45) that consists of zeros and ones) 
G=zeros(kbar,kbar*2);
CC=0;
for j=2:2:kbar*2;
    CC=CC+1;
   G(CC,(j-1):j)=1;
end

CC=0;D=zeros(kbar,kbar*2);
for j=2:2:(kbar*2-2);
    CC=CC+1;
D(CC,(j+1):kbar*2)=dst(1,1:kbar*2-j)*(gamma+varphi);
end;

D=D+G;

SigmaVV=zeros(kbar*4,kbar*4);
for j=1:2:kbar*2;
    SigmaVV(j,j)=sigma2v/(j^2);
    SigmaVV(j+1,j+1)=sigma2v/(j^2);
end;
SigmaVVst=SigmaVV;

W=zeros(kbar*4,kbar*4);
W(1:kbar*2,1:kbar*2)=Mst;
W(kbar*2+1:kbar*4,1:kbar*2)=eye(kbar*2);

Pst=SigmaVVst;
Mbar=zeros(kbar,kbar*2);
tic
% MAIN LOOP SOLUTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Diff=1;
while Diff > tol;
    
CC=0;
for j=1:2:kbar*2
    CC=CC+1;
Mbar(CC,(j+2:kbar*2))=cst(1,1:kbar*2-j-1)*M(1:kbar*2-j-1,1:kbar*2-j-1);
end;

c=a*D+b*Mbar;
cdiff=max(abs(c-cst));
cStep=Step*(1/(1+cdiff));
cst= cStep*c+ (1-cStep)*cst;    
    
d=dst*M-(fipi/gamma)*cst-(fiy/gamma)*dst+(1/gamma)*cst*Mst;
ddiff=max(abs(d-dst));
dst= (Step)*d + (1-Step)*dst;

CC=0;D=zeros(kbar,kbar*2);
for j=2:2:(kbar*2-2);
    CC=CC+1;
D(CC,(j+1):kbar*2)=dst(1,1:kbar*2-j)*(gamma+varphi);
end;
D=D+G;

L=zeros(3,4*kbar);
L(1,1:2)=1;
L(2,kbar*2+1:kbar*4)=cst;
L(3,kbar*2+1:kbar*4)=dst;


% %calculate P and K
Pdiff=1;

while Pdiff > tol;
    P=W*(Pst-Pst*L'*(inv(L*Pst*L'+SigmaEE))*L*Pst)*W'+SigmaVVst;
    Pdiff = max(  max( abs(Pst-P) )  );
    Pst=Step*P+(1-Step)*Pst;
end;

K = Pst*L'*(inv(L*Pst*L'+SigmaEE));

Nbar=K*L*N;
N(3:kbar*2,:)=Nbar(1:kbar*2-2,:)+[zeros(kbar*2-2,2),K(1:kbar*2-2,3);];

SigmaVV=N*SigmaUU*N';
VVdiff = max(  max( abs(SigmaVVst-SigmaVV) )  );
SigmaVVst=Step*SigmaVV+(1-Step)*SigmaVVst;
    

M(1,1)=rho;
Wbar=K*L*W;WWbar=W-K*L*W;
M(3:kbar*2,1:kbar*2)=[zeros(kbar*2-2,2),WWbar(1:kbar*2-2,1:kbar*2-2);]+[Wbar(1:kbar*2-2,1:kbar*2);];
Mdiff= max(  max( abs(Mst-M) )  );
Mst=Step*M+(1-Step)*Mst;
W(1:kbar*2,1:kbar*2)=Mst;
Diff=max([cdiff,ddiff,Mdiff,VVdiff]);
end;
%END OF MAIN SOLUTION LOOP%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
toc

%PLOTTING FIGURES ETC %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%Impulse response functions of inflationa and output%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
T=30;
shock=zeros(kbar,1);
for t=1:T;
    shock=M^(t-1)*N(1:kbar*2,1);         %use this line for persistent preference shock IRF
%     shock=M^(t-1)*N(1:kbar*2,2);       %use this line for transitory preference shock IRF
%     shock=M^(t-1)*N(1:kbar*2,3);       %use this line for demand shock IRF
    inf(t,1)=c*shock;
    output(t,1)=d*shock;
end;
hold on
figure(1)
plot(inf);
xlabel({'Inflation'},'FontSize',20);

hold on
figure(3)
plot(output);
xlabel({'Output'},'FontSize',20);

% Impulse reponse function marginal cost hierarchy%%%%%%%%%%%%%%%%%%%%%%%%%%
shock=zeros(kbar,1);
shock(1,1)=1;
for t=1:T;
    shock=D*M^(t-1)*N(1:kbar*2,1);           %use this line for persistent preference shock IRF
%     shock=M^(t-1)*N(1:kbar*2,2);           %use this line for transitory preference shock IRF
%     shock=M^(t-1)*N(1:kbar*2,3);           %use this line for output shock IRF
    mcimp(t,1)=shock(1,1);    
    mc1imp(t,1)=shock(2,1);
    mc2imp(t,1)=shock(3,1);
    mc3imp(t,1)=shock(4,1);    
    mc4imp(t,1)=shock(5,1);
    mc5imp(t,1)=shock(6,1);
 end;
figure(2);
hold on;
plot(mcimp,'k-');
plot(mc1imp,'k--');
plot(mc2imp,'k:');
plot(mc3imp,'b-');
xlabel({'Marginal cost hierarchy'},'FontSize',20);

%Simulate for Hybrid Phillips curve estimation%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
T=300;
qdata=zeros(2,T/3);
shockvar=[sigma2e;
      sigma2v;
      sigma2y;
      sigma2pi;];

for j=1:4;
shockserie(j,:)= randn(1,T)*(shockvar(j,1)^.5);
end;
piserie=zeros(1,T);
pl=zeros(1,T);
yserie=zeros(1,T);
mcserie=zeros(1,T);
H=zeros(kbar*2,T);
for i=1:T-1;
    H(:,i+1)=M*H(:,i)+N(1:kbar*2,:)*shockserie(2:4,i+1);
    piserie(1,i+1)=c*H(:,i+1);
    pl(1,i+1)=pl(1,i)+c*H(:,i+1);
    yserie(1,i+1)=d*H(:,i+1);
    mcserie(1,i+1)=[1,1]*H(1:2,i+1)+(gamma+varphi)*d*H(:,i+1);
end;

for i=1:3:T;
   p(1,(i+2)/3)=piserie(1,i)/3+piserie(1,i+1)/3+piserie(1,i+2)/3;
   y(1,(i+2)/3)=yserie(1,i)/3+yserie(1,i+1)/3+yserie(1,i+2)/3;
   mc(1,(i+2)/3)=mcserie(1,i)/3+mcserie(1,i+1)/3+mcserie(1,i+2)/3;
end;

qdata=[p;y;mc;]'; %Matrix with "quarterly" data (export to E-views or other GMM module)

%Simulate individual price path  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
calvo=rand(1,T);
shockser=[shockserie(1,:)+shockserie(2,:);
     shockserie(3,:);
      shockserie(4,:);];
  
pstarserie=zeros(1,T);
Hj=zeros(kbar*2,T);
absch=[];
for i=1:T-1;
    Hj(:,i+1)=M*Hj(:,i)+N(1:kbar*2,:)*shockser(:,i+1);
    if calvo(1,i+1) < theta
    pstarserie(1,i+1)=pstarserie(1,i);
    else
    pstarserie(1,i+1)=inv(1-theta)*c(1:kbar*2)*Hj(1:kbar*2,i+1)+(1-theta*beta)*shockserie(1,i+1)+pl(1,i);
    absch=[absch abs(pstarserie(1,i+1)-pstarserie(1,i))];           
    end
end

figure(4)
subplot(2,1,1);plot(piserie(1,1:100),'color','black','linewidth',2);
xlabel({'Inflation'},'FontSize',20);
subplot(2,1,2);plot(pl(1,1:100),'color','black','linewidth',2);
hold on
subplot(2,1,2);plot(pstarserie(1,1:100),'color','black','linewidth',2)
legend('Privcelevel','Price j');
xlabel({'Privcelevel and Price j'},'FontSize',20);

%compute relative average absolute price changes of price level and price of good j
display('Relative magnitude of indidual price changes over aggregate pricelevel changes ');
mean(absch)/mean(abs(piserie))




%compute variance ratio of idiosyncratic shock and mc
Hvar=dlyap(M,N(1:kbar*2,:)*SigmaUU*N(1:kbar*2,:)');
dd=(gamma+varphi)*d+[1,1,zeros(1,(kbar*2-2));];
mcvar=dd*Hvar*dd';
display('Idiosyncratic shock variance over average marginal cost variance ratio = ');
sigma2e/mcvar