% for the calibration, check gali p. 67-68
fi = 5;
ep = 9;
sig = 1.000000001;
alf = 0.25;
mu = ep/(ep-1);
rho_a = 0.85;
sig_a = 0.0095;
sig_v = 0.000625;
th = 3/4;
bet = 0.99;
taylor_pi = 1.5;
taylor_x = 0.5/4;
rho_v = 0.5;
rho_a = 0.9;

fi_ya = (1+fi)/(sig*(1-alf)+fi+alf);
%fi_y = -(1-alf)*(mu-log(1-alf))/(sig*(1-alf)+fi+alf);

lam = (1-th)*(1-bet*th)/th * (1-alf)/(1-alf+alf*ep);
kap = lam*(sig+(fi+alf)/(1-alf));

A = [rho_v          0                   0           0
        0           rho_a               0           0
        0           0               1/bet          -1/bet*kap
        1/bet    fi_ya*(1-rho_a)   (taylor_pi-1/bet)/sig     1+(taylor_x+kap/bet)/sig];
            
% solution to the model:  
[M,C] = ssolve(A,2); % 2 backward-looking variables, supply 2 as second argument
M = real(M); % remove tiny complex part due to roundig errors
C = real(C); 
            
            

e=zeros(2,16);   % order of shocks: a v 

% INSTRUCTION: uncomment ONE of the below lines to generate IRFS or a
% simulation
%e(1,1) = sig_v; % this is an mp shock
%e(1,1) = 0.25/100; % this is Galis MP shock

e(2,1) = sig_a; % this is a tech shock
% e = [sig_v 0; 0 sig_a] * randn(2,16);
n = 16; % desired lenght of sample/IRFS
x1 = zeros(2,n);
x2 = zeros(2,n);

for j=1:16
    x1(:,j) = e(:,j);
    if j>1
        x1(:,j) = x1(:,j)+M*x1(:,j-1);  % if j>1 add persistence component
    end
end

% Now we compute the other variables

x2 = C*x1;  % the forward-looking variables, ordered pi, x

y_n = fi_ya * x1(2,:); % natural output
r_n = -sig*fi_ya*(1-rho_a)*x1(2,:); % the natural real rate
    

tp = eye(2) - [ bet*rho_a           kap
                (rho_a  - taylor_pi)/sig/(1-taylor_x/sig)       rho_a/(1+taylor_x/sig)];
c_a = inv(tp)*[0;   -fi_ya*(1-rho_a)/(1+taylor_x/sig)];            

figure
subplot(2,3,1)
    plot(100*y_n);
    title('Natural output');
subplot(2,3,2)
    plot(400*r_n);
    title('Natural real rate');
    
subplot(2,3,3)
    plot(400*x2(1,:));
    title('Inflation, annualized, percentage points');
    
subplot(2,3,4)
    plot(100*x2(2,:));
    title('Output-gap, percentage points');

subplot(2,3,5)
    plot(100*x2(2,:)+100*y_n);
    title('Output, percentage points');
