% FULL INFORMATION PART OF BARSKY & SIMS MODEL

clear all
close all

%%%%
beta=0.99;
eta=1;
alpha=0.33;
delta=0.025;
s=1;
rhoa=1;
rhog=0.95;
mu=1.1;
theta=0.55;
gamma=0.15;
rhoi=0.75;
phi1=0.5;
phi2=0.36;
rhoga=0.75;
psi=0;
rhoz=0.9;
vv=0.66;
gstar=0.2;

nimpstep = 1000;	          
numplot = 30;
per=100;                    %%number of periods per year

% Calibration

rstar = 1/beta - (1 - delta); % steady state interest rate
kstar = (alpha/(rstar*mu))^(1/(1 - alpha)); % steady state capital per worker
%gstar = 0.2; % steady state government consumption percentage
Nstar = ((1/mu)*((1-alpha)*kstar^(alpha))*((1-gstar)*kstar^(alpha)-delta*kstar)^(-1))^(1/(1+1/eta)); % steady state employment
cstar = Nstar*(kstar^(alpha)-delta*kstar-gstar*kstar^(alpha)); % steady state consumption
Kstar = Nstar*kstar; % steady state capital
Ystar = Kstar^(alpha)*Nstar^(1-alpha); % steady state output
Gstar = gstar*Ystar; % steady state government consumption level
Istar = Ystar - cstar - Gstar; % steady state investment
Istar2 = delta*Kstar; % check
v = (1-delta)/(alpha*(1/mu)*(Ystar/Kstar) + (1-delta));

    

%******************************************************************
%**** Begin setup for constructing log-linear system coef matrix **
%******************************************************************

nlead = 1;  % Number of leads in system 
nlag  = 1;  % Number of lags in system 

% Coefficient indicators: enter the name of the variables, including the innovation to exogenous variable
%(note xnames is not used for anything right now except to determine 
% number of equations "neq" )
      
               % 
               
xnames = strvcat('c','a','k','N','I','g','y','G','w','r','q','i','pi','p','mc','cf','kf','yf','Nf','If','wf','rf','qf','e5y','ga','z','lam','ea','eg','ega','ei','ee','ez');


xnum = length(xnames); %  Number of variables in system 
neq = xnum;            %  Number of equations (same)

nlag1=nlag;
nlead1=nlead;


% Ordering variables in the X vector using position indicators 
                                                        
% Note: All shock variables must be placed last in the system!

cpos = 1;
apos = 2;
kpos = 3;
Npos = 4;
Ipos = 5;
gpos = 6;
ypos = 7;
wpos = 8;
Gpos = 9;
rpos = 10;
qpos = 11;
ipos = 12;
pipos = 13;
ppos = 14;
mcpos = 15;
cfpos = 16;
kfpos = 17;
yfpos = 18;
Nfpos = 19;
Ifpos = 20;
wfpos = 21;
rfpos = 22;
qfpos = 23;
e5ypos = 24;
gapos = 25;
zpos = 26;
eapos = 27;
egpos = 28;
egapos = 29;
eipos = 30;
lampos = 31;
eepos = 32;
ezpos = 33;
%ecpos = 29;


nlag2=nlag;
nlead2=nlead;


% these are position indicators for leads and lags in constructing
% the coefficient matrix 

colzero = 0+nlag*xnum;      % Position counter for start of contemp. coefs 
			                % i.e. if there is one lag and 11 equations (including the innovation) then 
                            % the "first" contemporaneous variable is at
                            % position 12
collead = 0+nlag*xnum+xnum; % Position counter for start of lead coefs 
collag  = 0;                % Position counter for start of lag coefs  

% Indicators for contemporanous coefs for each variable: 

nlag3=nlag;
nlead3=nlead;


azero = colzero + apos;
yzero = colzero + ypos;
Nzero = colzero + Npos;
Izero = colzero + Ipos;
kzero = colzero + kpos;
czero = colzero + cpos;
gzero = colzero + gpos;
Gzero = colzero + Gpos;
wzero = colzero + wpos;
rzero = colzero + rpos;
qzero = colzero + qpos;
izero = colzero + ipos;
pizero = colzero + pipos;
pzero = colzero + ppos;
mczero = colzero + mcpos;
cfzero = colzero + cfpos;
kfzero = colzero + kfpos;
yfzero = colzero + yfpos;
Nfzero = colzero + Nfpos;
Ifzero = colzero+ Ifpos;
wfzero = colzero + wfpos;
rfzero = colzero + rfpos;
qfzero = colzero + qfpos;
zzero = colzero + zpos;
gazero = colzero + gapos;
e5yzero = colzero + e5ypos;
eazero = colzero + eapos;
egzero = colzero + egpos;
egazero = colzero + egapos;
eizero = colzero + eipos;
eezero = colzero + eepos;
ezzero = colzero + ezpos;
lamzero = colzero + lampos;
%eczero = colzero + ecpos;



% Indicators for lead coefficients for each variable: 
nlag4=nlag;
nlead4=nlead;

alead = collead + apos;
ylead = collead + ypos;
Nlead = collead + Npos;
Ilead = collead + Ipos;
klead = collead + kpos;
clead = collead + cpos;
glead = collead + gpos;
Glead = collead + Gpos;
wlead = collead + wpos;
rlead = collead + rpos;
qlead = collead + qpos;
ilead = collead + ipos;
pilead = collead + pipos;
plead = collead + ppos;
mclead = collead + mcpos;

cflead = collead + cfpos;
kflead = collead + kfpos;
yflead = collead + yfpos;
Nflead = collead + Nfpos;
Iflead = collead+ Ifpos;
wflead = collead + wfpos;
rflead = collead + rfpos;
qflead = collead + qfpos;
galead = collead + gpos;
e5ylead = collead + e5ypos;
lamlead = collead + lampos;
ealead = collead + eapos;
eglead = collead + egpos;
egalead = collead + egapos;
eilead = collead + eipos;
eelead = collead + eepos;
ezlead = collead + ezpos;
zlead = collead + zpos;
%eclead = collead + ecpos;



% Indicators for lag coefficients for each variable: 
nlag5=nlag;
nlead5=nlead;

alag = collag + apos;
ylag = collag + ypos;
Nlag = collag + Npos;
Ilag = collag + Ipos;
klag = collag + kpos;
clag = collag + cpos;
glag = collag + gpos;
Glag = collag + Gpos;
wlag = collag + wpos;
qlag = collag + qpos;
rlag = collag + rpos;
ilag = collag + ipos;
pilag = collag + pipos;
plag = collag + ppos;
mclag = collag + mcpos;
e5ylag = collag + e5ypos;
cflag = collag + cfpos;
kflag = collag + kfpos;
yflag = collag + yfpos;
Nflag = collag + Nfpos;
Iflag = collag + Ifpos;
wflag = collag + wfpos;
rflag = collag + rfpos;
qflag = collag + qfpos;
galag = collag + gapos;
zlag = collag + zpos;
ealag = collag + eapos;
eglag = collag + egpos;
egalag = collag + egapos;
eilag = collag + eipos;
eelag = collag + eepos;
ezlag = collag + ezpos;
lamlag = collag + lampos;
%eclag = collag + ecpos;


nlag6=nlag;
nlead6=nlead;


% Now we have one vector with all of the leads, contempraneous and lags stacked in one column
% Determine number of total posible coefficients per equation: 

ncoef = neq*(nlag+nlead+1);

%Initialize the coeffienct matrix, where each row is an equation of the model

cof = zeros(neq,ncoef);      % Coef matrix --- Each row is an equation



% Setup coefficients vectors for each equation: 
% ==============================================

% (1) Euler equation: c(t+1) - c(t) - s*(r(t)) = 0;
cof(1,lamlead) = 1;
cof(1,lamzero) = -1;
cof(1,rzero) = 1;
%cof(1,zzero) = 1;

% (2) Accounting condition: y(t) = (c/y)*c(t) + (i/y)*i(t) + (g/y)*G(t)
cof(2,yzero) = 1;
cof(2,czero) = -(cstar/Ystar);
cof(2,Izero) = -(Istar/Ystar);
cof(2,Gzero) = -gstar;

% (3) Production: y(t) = a(t) + alpha*k(t) + (1-alpha)*n(t)
cof(3,yzero) = 1;
cof(3,azero) = -1;
cof(3,Nzero) = -(1-alpha);
cof(3,klag) = -alpha;

% (4) Labor demand: w(t) = mc(t) + a(t) + alpha*k(t) - alpha*N(t)
cof(4,wzero) = 1;
cof(4,azero) = -1;
cof(4,klag) = -alpha;
cof(4,Nzero) = alpha;
cof(4,mczero) = -1;

% (5) Capital demand (1-v)*(mc(t+1) + y(t+1) - k(t)) + v*g(t+1) - q(t) =
% r(t)
cof(5,rzero) = 1;
cof(5,mclead) = -(1-v);
cof(5,ylead) = -(1-v);
cof(5,kzero) = (1-v);
cof(5,qlead) = -v;
cof(5,qzero) = 1;

% (5) Capital demand: R(t) = a(t) + (alpha-1)*k(t) + (1-alpha)*N(t)
%cof(5,Rzero) = 1;
%cof(5,azero) = -1;
%cof(5,klag) = -(alpha-1);
%cof(5,Nzero) = -(1-alpha);

% (6) Fisher relationship: r(t) = i(t) - pi(t+1)
cof(6,rzero) = -1;
%cof(6,rlead) = -1;
cof(6,izero) = 1;
cof(6,pilead) = -1;

% (7) Labor supply: (1/eta)*n(t) = -(1/s)*c(t) + w(t)
cof(7,Nzero) = (1-psi)*1/eta;
cof(7,lamzero) = -(1-psi)*1/s;
cof(7,wlag) = psi;
cof(7,wzero) = -1;

% (8) Process for technology: a(t) = rhoa*a(t-1) + ea(t)
cof(8,azero) = 1;
cof(8,galag) = -1;
cof(8,alag) = -rhoa;
cof(8,eazero) = -1;

%(9) Relationship between g and G: G(t) = g(t) + y(t)
cof(9,Gzero) = 1;
cof(9,gzero) = -1;
cof(9,yzero) = -1;

% (10) Process for g: g(t) = rhog*g(t-1) + eg(t)
cof(10,gzero) = 1;
cof(10,glag) = -rhog;
cof(10,egzero) = -1;

% (11) Accumulation equation k(t) = delta(i) + (1-delta)*k(t-1)
cof(11,kzero) = 1;
cof(11,Izero) = -delta;
cof(11,klag) = -(1-delta);

% (12) Definition of q: q(t) = gamma*(i(t) - k(t-1))
cof(12,qzero) = -1;
cof(12,Izero) = gamma;
cof(12,klag) = -gamma;

% (13) Phillips Curve: pi(t) = mc(t)*(1-theta)*(1-theta*beta)/theta +
% beta*pi(t+1)
cof(13,pizero) = -1;
cof(13,pilead) = beta;
cof(13,mczero) = (1-theta)*(1-theta*beta)/theta;
%cof(13,eczero) = 1;

% (14) Definition of inflation and the price level: pi(t) = p(t) - p(t-1)
cof(14,pizero) = -1;
cof(14,pzero) = 1;
cof(14,plag) = -1;

% (15) Monetary policy rule: i(t) = rhoi*i(t-1) + phi1*pi(t)
cof(15,izero) = -1;
cof(15,ilag) = rhoi;
cof(15,pizero) = phi1;
%cof(15,pilead) = phi1;
cof(15,yzero) = phi2;
%cof(15,yfzero) = -phi2;
cof(15,ylag) = -phi2;
cof(15,eizero) = 1;
%cof(15,ylag) = -phi2;


% Now do flexible price versions

cof(16,cflead) = 1;
cof(16,cfzero) = -1;
cof(16,rfzero) = -s;
%cof(1,rlead) = -s;

% (2) Accounting condition: y(t) = (c/y)*c(t) + (i/y)*i(t) + (g/y)*G(t)
cof(17,yfzero) = 1;
cof(17,cfzero) = -(cstar/Ystar);
cof(17,Ifzero) = -(Istar/Ystar);
cof(17,Gzero) = -gstar;

% (3) Production: y(t) = a(t) + alpha*k(t) + (1-alpha)*n(t)
cof(18,yfzero) = 1;
cof(18,azero) = -1;
cof(18,Nfzero) = -(1-alpha);
cof(18,kflag) = -alpha;

% (4) Labor demand: w(t) = mc(t) + a(t) + alpha*k(t) - alpha*N(t)
cof(19,wfzero) = 1;
cof(19,azero) = -1;
cof(19,kflag) = -alpha;
cof(19,Nfzero) = alpha;

% (5) Capital demand (1-v)*(mc(t+1) + y(t+1) - k(t)) + v*g(t+1) - q(t) =
% r(t)
cof(20,rfzero) = 1;
cof(20,yflead) = -(1-v);
cof(20,kfzero) = (1-v);
cof(20,qflead) = -v;
cof(20,qfzero) = 1;

% (11) Accumulation equation k(t) = delta(i) + (1-delta)*k(t-1)
cof(21,kfzero) = 1;
cof(21,Ifzero) = -delta;
cof(21,kflag) = -(1-delta);

% (12) Definition of q: q(t) = gamma*(i(t) - k(t-1))
cof(22,qfzero) = -1;
cof(22,Ifzero) = gamma;
cof(22,kflag) = -gamma;

% (7) Labor supply: (1/eta)*n(t) = -(1/s)*c(t) + w(t)
cof(23,Nfzero) = 1/eta;
cof(23,cfzero) = 1/s;
cof(23,wfzero) = -1;

cof(24,gazero) = -1;
cof(24,galag) = rhoga;
cof(24,egazero) = 1;

cof(25,e5yzero) = -1;
cof(25,e5ylag) = 0.8;
cof(25,egzero) = 0.05;
cof(25,eazero) = 0.5;
cof(25,egazero) = 4.5;
cof(25,eezero) = 1;

cof(26,zzero) = -1;
cof(26,zlag) = rhoz;
cof(26,ezzero) = 1;

cof(27,lamzero) = -1;
cof(27,czero) = -((1/((1-vv)*(1-vv*beta)))*(1 + beta*vv^2));
cof(27,clag) = (1/((1-vv)*(1-vv*beta)))*vv;
cof(27,clead) = (1/((1-vv)*(1-vv*beta)))*vv*beta; % not ylead, perceived tomorrow's y!
cof(27,zzero) = 1/(1-beta*vv);
cof(27,zlead) = -beta*vv/(1-beta*vv);

% (12) Shock process for technology
cof(28,eazero) = 1;

% (13) Shock for government spending
cof(29,egzero) = 1;

cof(30,egazero) = 1;

cof(31,eizero) = 1;

cof(32,eezero) = 1;

cof(33,ezzero) = 1;

%cof(29,eczero) = 1;

[DD,dd] = size(cof);

Q = isnan(cof);

Q3 = max(cof);

for j = 1:dd
    QQ(:,j) = max(Q(:,j));
end

%if max(QQ)<1&&max(Q3)<9999

%*****************************************************************
%*************Begin Solution Algorithm ***************************
%*****************************************************************
%                                                                 
%  Solve a linear perfect foresight model using the gauss eig     
%  function to find the invariant subspace associated with the big 
%  roots.  This procedure will fail if the companion matrix is     
%  defective and does not have a linearly independent set of       
%  eigenvectors associated with the big roots.                     
%                                                                  
%  Input arguments:                                                
%                                                                  
%    h         Structural coef matrix (neq,neq*(nlead+1+nlag)).     
%    neq       Number of equations.                                
%    nlead     Number of leads.                                     
%    nlag      Number of lags.                                      
%    condn     lag tolerance used as a condition number test       
%              by numeric_shift and reduced_form.                  
%    uprbnd    Inclusive upper bound for the modulus of roots      
%              allowed in the reduced form.                        
%                                                                  
%  Output arguments:                                               
%                                                                  
%    cofb      Reduced form coefficient matrix (neq,neq*nlag).     
%    scof      Observable Structure                                
%    amat      Companion form matrix                               
%    b         Contemporaneous coefficient matrix                  

%    Model satisfies:                                              
%                                                                  
%    z(t) = amat*z(t-1) + b*e(t)                                   
%                                                                  
%    where the first neq elements of z(t) are the contemporaneous  
%    values of the variables in the model                          
%    and e(t) is the shock vector of conformable dimension         
%                                                                  
%    rts       Roots returned by eig.                              
%    ia        Dimension of companion matrix (number of non-trivial
%              elements in rts).                                   
%    nexact    Number of exact shiftrights.                        
%    nnumeric  Number of numeric shiftrights.                      
%    lgroots   Number of roots greater in modulus than uprbnd.     
%    mcode     Return code: see function aimerr.                   
%*****************************************************************


% Use AIM procedure to solve model: 
uprbnd = 1+1e-8;    % Tolerance values for AIM program 
condn = 1e-8;

% ---------------------------------------------------------------------
% Run AIM
% ---------------------------------------------------------------------

[cofb,rts,ia,nex,nnum,lgrts,mcode] = ...
       SPAmalg(cof,neq,nlag,nlead,condn,uprbnd);

scof = SPObstruct(cof,cofb,neq,nlag,nlead);

% need to calculate amat and b
% ===============================
s0 = scof(:,(neq*nlag+1):neq*(nlag+1)); %Contemp. coefs from obs.
                                        %structure
amat=zeros(neq*nlag,neq*nlag);   		% Initialize A matrix 
bmat=cofb(1:neq,((nlag-1)*neq+1):nlag*neq); % Lag 1 coefficients 
i=2;
while i<=nlag;
  bmat=[bmat cofb(1:neq,((nlag-i)*neq+1):(nlag-i+1)*neq)]; % Lag i coefs 
  i=i+1;
end;
amat(1:neq,:)=bmat;  				% Coefs for equations 
if nlag>1;
 amat((length(cofb(:,1))+1):length(amat(:,1)),1:neq*(nlag-1))=eye(neq*(nlag-1));
end;
b = zeros(length(amat(:,1)),length(s0(1,:)));
b(1:length(s0(:,1)),1:length(s0(1,:))) = inv(s0);  % Store coefs