/* **************************************** */
/* Gauss Code to Generate IRF for RBC Model */
/* **************************************** */

/* Parameter Values */
theta = 0.64;
delta = 0.025;
phi = 2/3;
rho = 0.95;
beta = 0.99;
sigma = 0.763;
nobs = 50;

/* Steady-State Ratios */
yratk = (1/(1-theta))*((1/beta)-(1-delta));
kraty = 1/yratk;
cratk = yratk - delta;
craty = cratk/yratk;
nbar = 1/(1+(phi*craty/(theta*(1-phi))));

/* Intermediate Parameter Values */
denom = (1-nbar)*theta - 1;
a1 = (1-nbar)*theta/denom;
a2 = -(1-theta)/denom;
a3 = -1/denom;
a4 = a1 - craty;
a5 = a2 + kraty*(1-delta);
a6 = 1 - (beta*(1-theta)*yratk*a1);
a7 = beta*(1-theta)*yratk*(1-a2);
a8 = -beta*(1-theta)*yratk*a3;

/* Define Matrices */
F = zeros(3,3);
F[1,1] = a5;
F[1,2] = a4;
F[1,3] = a3;
F[2,2] = 1;
F[3,3] = rho;

G = zeros(3,3);
G[1,1] = kraty;
G[2,1] = a7;
G[2,2] = a6;
G[2,3] = a8;
G[3,3] = 1;

H = zeros(3,4);
H[2,2] = a7;
H[2,3] = a6;
H[2,4] = a8;
H[3,1] = -1;

/* Calculate Eigenvalues and Eigenvectors */
A = inv(F)*G;
{eigs,eigvecs} = eigv(A);
qinv = inv(eigvecs);

/* VAR Coefficient Matrix */
b1 = yratk*(a5 - a4*qinv[1,1]/qinv[1,2]);
b2 = yratk*(a3 - a4*qinv[1,3]/qinv[1,2]);
ar1coef = zeros(2,2);
ar1coef[1,1] = b1;
ar1coef[1,2] = b2;
ar1coef[2,2] = rho;
print "State VAR(1) Matrix = " ar1coef;

/* Policy Function Matrices */
pimat1 = zeros(3,3);
pimat1[1,1] = 1;
pimat1[1,2] = -a1;
pimat1[2,2] = qinv[1,2];
pimat1[3,1] = 1;
pimat1[3,2] = -1;
pimat1[3,3] = -1/(1-nbar); 

pimat2 = zeros(3,2);
pimat2[1,1] = a2;
pimat2[1,2] = a3;
pimat2[2,1] = -qinv[1,1];
pimat2[2,2] = -qinv[1,3];
pimat = inv(pimat1)*pimat2;
print "Pi Matrix = " pimat;

/* ********************************* */
/* Calculating and Graphing the IRFs */
/* ********************************* */

eps = zeros(nobs,1);
k = zeros(nobs,1);
z = zeros(nobs,1);
y = zeros(nobs,1);
c = zeros(nobs,1);
n = zeros(nobs,1);
eps[1] = sigma;

z[1] = eps[1];
for ii(2,nobs,1);
    z[ii] = rho*z[ii-1] + eps[ii];
   	k[ii] = b1*k[ii-1] + b2*z[ii-1];
endfor;

for ii(1,nobs,1);
	y[ii] = pimat[1,1]*k[ii] + pimat[1,2]*z[ii];
	c[ii] = pimat[2,1]*k[ii] + pimat[2,2]*z[ii];
	n[ii] = pimat[3,1]*k[ii] + pimat[3,2]*z[ii];
endfor;

x = seqa(1,nobs,1);
irfs = k~y~c~n;
library pgraph;
graphset;
title("Impulse Response Function for a Simple RBC Model");
_plegctl = 1;
_plegstr = "K\000Y\000C\000N";
xlabel("Time");
xy(x,irfs);