/* **************************************** */
/* 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;

/* 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 for Procedure */
A = zeros(3,6);
A[1,1] = a5;
A[1,4] = -kraty;
A[1,5] = a4;
A[1,6] = a3;
A[2,4] = -a7;
A[2,5] = 1;
A[3,3] = rho;
A[3,6] = -1;

B = zeros(3,3);
B[2,2] = -a6;
B[2,3] = -a8; 

C = zeros(3,1);
C[3,1] = 1;

D = 0;
F = 0;

G = zeros(2,6);
G[1,1] = a2;
G[1,5] = a1;
G[1,6] = a3;
G[2,1] = (1-nbar)*a2;
G[2,5] = (1-nbar)*(a1-1);
G[2,6] = (1-nbar)*a3;

/* Parameters for Procedure */
sig = {0.763};
K = 75;
L = 0;
T = 75;

/* Run Procedure */
{irf,std,cor} = resolve(A,B,C,D,F,G,sig,K,L,T);