%%%
%%% Long-Run Risk and News Risk in a RBC model
%%%

clear all ; close all ; clc ;

%%%
%%% parameters
%%%

gamma = 1 ; %%% 1/IES = RA
epshw = 1e-6 ; %%% Frisch elasticity (if close to zero, this is the fixed labor case)
nstar = 0.2 ; %%% SS employment
delta = 0.08 ; %%% depreciation rate
mu = 0.02 ; %%% trend TFP growth
alpha = .3 ; %%% capital share
eta = 1e-6; %%% adj cost parameter (elasticity of I/K to Tobin's q) -> close to 0, no adjustment cost
beta = .99 ; %%% discount factor

showall = 0 ; %%% if set to 1, this shows some additional graphs

%%%
%%% shock process to choose
%%% note: for now, cannot choose both rhoz>0 and news:
%%% it's either one or the other.

lagnews = 0 ;
%%% number of lags. (note: it should be easy to change the program to have dlogZ(t) depending on various lags,
%%% e.g. eps(t-4), eps(t-5) ,etc., a smooth diffusion)
rhoz = 0 ; %%% persistence of growth rate 

%%%
%%% steady-state
%%%

kovern = ((1/beta*exp(mu*gamma)-1+delta)/alpha)^(1/(alpha-1))*exp(mu) ;
kstar = kovern*nstar ;
ystar = (kstar*exp(-mu))^alpha*nstar^(1-alpha) ;
istar = (1-(1-delta)*exp(-mu))*kstar ;
cstar = ystar - istar ;
qstar = 1 ;

%%%
%%% elasticities
%%%

epshlambda = 1/gamma*epshw ;
epscw = epshw*(gamma-1)/gamma*(1-alpha)/(cstar/ystar) ;
epsclambda = (epscw-1)/gamma ;

%%%
%%% define variables
%%%

currentindex = 1 ;
indexK = currentindex ; %%% this is K_t / z_t-1
currentindex = currentindex + 1 ;

if lagnews==0
    indexDLOGZ = currentindex ;
    currentindex = currentindex + 1 ;
end

if lagnews>0
    for i=1:lagnews
        aa = strcat('indexEPS',num2str(i),'=','currentindex;');
        eval(aa);
        currentindex = currentindex + 1 ;
    end
    indexDLOGZ = currentindex ; %%% not a state any more.
    currentindex = currentindex + 1 ;
end

if rhoz~=0;
    indexDLOGZlag = currentindex ;
    currentindex = currentindex + 1 ;
end

indexC = currentindex ;
currentindex = currentindex + 1 ;

indexLambda = currentindex ;
currentindex = currentindex + 1 ;

indexN = currentindex;
currentindex = currentindex + 1 ;

indexI = currentindex;
currentindex = currentindex + 1 ;

indexY = currentindex;
currentindex = currentindex + 1 ;

indexW = currentindex;
currentindex = currentindex + 1 ;

indexQ = currentindex;

if rhoz==0
    listvar = {'K/z(-1)','dlogZ','C/z','Lamba*Z(gamma)','N','I/z','Y/z','W/z','Tobin Q'};
end
if rhoz>0
    listvar = {'K/z(-1)','dlogZ','dlogZ(-1)','C/z','Lamba*Z(gamma)','N','I/z','Y/z','W/z','Tobin Q'};
end
if lagnews>0
    listvar = {'K/z(-1)'} ;
    for i=1:lagnews;
        listvar{1+i} = strcat('indexEPS',num2str(i)) ;
    end
    listvar{1+lagnews+1} = 'dlogZ';
    listvar{2+lagnews+1} = 'C/z';
    listvar{2+lagnews+1+1} = 'Lamba*Z(gamma)';
    listvar{2+lagnews+2+1} = 'N' ;
    listvar{2+lagnews+3+1} = 'I/z' ;
    listvar{2+lagnews+4+1} = 'Y/z' ;
    listvar{2+lagnews+5+1} = 'W/z';
    listvar{2+lagnews+6+1} = 'Tobin Q';
end

%%%
%%% now stack the log-linearized system:
%%%

Arbc = zeros(currentindex);
Brbc = zeros(currentindex);


linen = 1 ;

%%% "shocks" for Z
if lagnews==0

    if rhoz==0
        Arbc(linen,indexDLOGZ) = 1;
    end

    if rhoz~=0
        Arbc(linen,indexDLOGZ) = 1;
        Brbc(linen,indexDLOGZ)= rhoz ;

        linen = linen + 1 ;
        Arbc(linen,indexDLOGZlag) = 1 ;
        Brbc(linen,indexDLOGZ) = 1 ;
    end

end

if lagnews>0
    
    Brbc(linen, indexDLOGZ) = -1;
    Brbc(linen, indexEPS1) = 1;
    
    for i=1:lagnews-1
        linen = linen + 1 ;
        uu = strcat('indexEPS',num2str(i)) ;
        vv = strcat('indexEPS',num2str(i+1)) ;
        eval(strcat('Arbc(linen,',uu,')=1'));
        eval(strcat('Brbc(linen,',vv,')=1'));
    end
    
    i = lagnews ;
    linen = linen + 1 ;
    uu = strcat('indexEPS',num2str(i)) ;
    eval(strcat('Arbc(linen,',uu,')=1'));

end

%%% Express C in function of Lambda and W using epshw, gamma:
linen = linen + 1 ;
Brbc(linen,indexLambda) = epsclambda ;
Brbc(linen,indexW) = epscw ;
Brbc(linen,indexC) = - 1 ;

%%% Express N in function of Lambda and W using epshw, gamma:
linen = linen + 1 ;
Brbc(linen,indexN) = -1 ;
Brbc(linen,indexW) = epshw ;
Brbc(linen,indexLambda) = epshlambda;

%%% Q-theory:
linen = linen + 1 ;
Brbc(linen, indexQ) = -1 ;
Brbc(linen, indexI) = eta ;
Brbc(linen, indexK) = -eta ;
Brbc(linen, indexDLOGZ) = eta ;

%%% Resource Constraint:
linen = linen + 1 ;
Brbc(linen, indexC) = cstar/ystar ;
Brbc(linen, indexI) = istar/ystar ;
Brbc(linen, indexY) = -1; %%% this is Y/z

%%% Production Function:
linen = linen + 1 ;
Brbc(linen, indexY) = -1 ;
Brbc(linen, indexK) = alpha;
Brbc(linen, indexDLOGZ) = -alpha;
Brbc(linen, indexN) = 1-alpha;

%%% Capital accumulation:
linen = linen + 1 ;
Arbc(linen, indexK) = 1 ;
Brbc(linen, indexK) = (1-delta)*exp(-mu) ;
Brbc(linen, indexI) = 1 - (1-delta)/exp(mu) ;
Brbc(linen, indexDLOGZ) = - (1-delta)/exp(mu) ;

%%% Euler equation:
linen = linen + 1 ;
Brbc(linen, indexLambda) = 1;
Brbc(linen, indexQ) = 1;
Arbc(linen, indexLambda) = 1;
Arbc(linen, indexQ) = (1-delta)*qstar/( (1-delta)*qstar + alpha*(kovern/exp(mu))^(alpha-1) ) ;
Arbc(linen, indexDLOGZ) = - gamma;
Arbc(linen, indexK) = (alpha*(alpha-1)*(kstar/nstar/exp(mu))^(alpha-1)) /( (1-delta)*qstar + alpha*(kstar/nstar/exp(mu))^(alpha-1) ) ;
Arbc(linen, indexN) = - Arbc(linen, indexK) ;
Arbc(linen, indexDLOGZ) = Arbc(linen, indexDLOGZ) - Arbc(linen, indexN);

%%% Def of wage as MPL
linen = linen + 1;
Brbc(linen, indexW) = -1 ;
Brbc(linen, indexK) = alpha ;
Brbc(linen, indexN) = -alpha ;
Brbc(linen, indexDLOGZ) = -alpha ;

if lagnews==0
    if rhoz==0
        numState = 2 ;

    else
        numState = 3 ;
    end
end

if lagnews>0
    numState = 1+lagnews ;
end
[M,N] = solab(Arbc,Brbc,numState);

%%%
%%% plot IRF
%%%

listinterestIR= (1:currentindex) ; %%% IRFs to plot
T = 30 ; %%% length of plot
if rhoz==0
    impulse = [0,1] ;
end
if rhoz>0
    impulse = [0,1,0] ;
end
if lagnews>0
    impulse = [0,zeros(1,lagnews-1),1];
end

imprep1=zeros(length(listinterestIR),T);
for k=1:T
    imprep1(1:numState,k)=N^(k-1)*impulse';
    imprep1(numState+1:end,k)=M*N^(k-1)*impulse' ;
end

if showall==1
    for i=1:length(listinterestIR);
        figure;
        plot(imprep1(listinterestIR(i),:),'b+-');
        title(['IR of :  ',char(listvar(listinterestIR(i)))]);
        grid on;
    end
end

%%% Now plot the 'real' variables'
Z = cumsum(imprep1(indexDLOGZ,:));
C = imprep1(indexC,:) + Z ;
I = imprep1(indexI,:) + Z ;
Y = imprep1(indexY,:) + Z ;
N = imprep1(indexN,:) ;
K = imprep1(indexK,2:end)+ Z(1:end-1) ;
K = [0,K];
Q = imprep1(indexQ,:);
W = imprep1(indexW,:) + Z ;

figure; plot((0:T-1),Z,'bo-','LineWidth',2)
title('IRF of Z');xlabel('periods after shock');ylabel('% deviation from steady-state'); grid on;
figure; plot((0:T-1),imprep1(indexDLOGZ,:),'bo-','LineWidth',2)
title('IRF of DLOGZ');xlabel('periods after shock');ylabel('% deviation from steady-state'); grid on;
figure; plot((0:T-1),Y,'bo-','LineWidth',2)
title('IRF of Y');xlabel('periods after shock');ylabel('% deviation from steady-state'); grid on;
figure; plot((0:T-1),N,'bo-','LineWidth',2)
title('IRF of N');xlabel('periods after shock');ylabel('% deviation from steady-state'); grid on;
figure; plot((0:T-1),C,'bo-','LineWidth',2)
title('IRF of C');xlabel('periods after shock');ylabel('% deviation from steady-state'); grid on;
figure; plot((0:T-1),I,'bo-','LineWidth',2)
title('IRF of I');xlabel('periods after shock');ylabel('% deviation from steady-state'); grid on;
figure; plot((0:T-1),K,'bo-','LineWidth',2)
title('IRF of K');xlabel('periods after shock');ylabel('% deviation from steady-state'); grid on;
figure; plot((0:T-1),Q,'bo-','LineWidth',2)
title('IRF of Q');xlabel('periods after shock');ylabel('% deviation from steady-state'); grid on;
figure; plot((0:T-1),W,'bo-','LineWidth',2)
title('IRF of W');xlabel('periods after shock');ylabel('% deviation from steady-state'); grid on;

%%%
%%% can plot things differently, e.g. several series on the same graph...
%%%