%RBC_solve_sim_simpl
%Solves and simulates RBC model with government spending
%Simplified version
%Uses function RE_LIN_SOLVE.M
%Economics 210C, Spring 2019, Garey Ramey, UCSD

%SET PARAMETER VALUES 
alpha = .34;
beta = 1.04^(-.25);
rho = .95;
sigma_eps = .007;
chi = 1;
eta = 10;
delta = .01;
G = .9;
gamma = .95;
sigma_nu = .025;
        
%CALCULATE STEADY STATE
kappa = ((1/beta-1+delta)/alpha)^(1/(alpha-1));
L_fun = @(x) (1-alpha)*kappa^alpha*x^(-1/eta)/chi ...
    +(delta*kappa-kappa^alpha)*x+G;
L = fsolve(@(x) L_fun(x),.2,optimoptions('fsolve','Display','off'));
C = (1-alpha)*kappa^alpha*L^(-1/eta)/chi;
K = kappa*L;
Y = kappa^alpha*L;
ratios = [C/Y G/Y];       %Steady state ratios

%CALCULATE RE SOLUTION
%Set up matrices of vector system
Av = [1 0; 0 -(alpha+1/eta)];
Alam = [1; 1];
Ak = [0; alpha];
Az = [0 0; 1 0];
Bv = [0 beta*alpha*(1-alpha)*kappa^(alpha-1)];
Blam1 = 1; 
Bk = -beta*alpha*(1-alpha)*kappa^(alpha-1); 
Bz = [beta*alpha*kappa^(alpha-1) 0];
Blam0 = 1;
Ck1 = kappa*L; 
Cv = [-C (1-alpha)*kappa^alpha*L]; 
Ck0 = (alpha*kappa^(alpha-1)+1-delta)*kappa*L; 
Cz = [kappa^alpha*L -G];
PI = [rho 0; 0 gamma];
    
Amat = {Av Alam Ak Az};
Bmat = {Bv Blam1 Bk Bz Blam0};
Cmat = {Ck1 Cv Ck0 Cz};

%Solve for REE using Sims' method
[PHI_U,PHI_S,B_S,PHI_v] = REE_lin_solve(Amat,Bmat,Cmat,PI);

%GENERATE ONE REPLICATION OF SIMULATED DATA
%Initialize
T_init = 200;  
T_samp = 4000;
T_sim = T_init+T_samp;
k_sim = zeros(T_sim,1);
z_sim = zeros(T_sim,1);
g_sim = zeros(T_sim,1);
lam_sim = zeros(T_sim,1);
c_sim = zeros(T_sim,1);
l_sim = zeros(T_sim,1);

%Draw shocks, calculate solution path
shk = [randn(T_sim,1)*sigma_eps randn(T_sim,1)*sigma_nu];
for t = 2:T_sim
    k_sim(t) = PHI_S(1,:)*[k_sim(t-1) z_sim(t-1) g_sim(t-1)]';
    z_sim(t) = PHI_S(2,:)*[k_sim(t-1) z_sim(t-1) g_sim(t-1)]'+shk(t,1);
    g_sim(t) = PHI_S(3,:)*[k_sim(t-1) z_sim(t-1) g_sim(t-1)]'+shk(t,2);
    lam_sim(t) = PHI_U*[k_sim(t) z_sim(t) g_sim(t)]';
    c_sim(t) = PHI_v(1,:)*[lam_sim(t) k_sim(t) z_sim(t) g_sim(t)]';
    l_sim(t) = PHI_v(2,:)*[lam_sim(t) k_sim(t) z_sim(t) g_sim(t)]';
end

%Calculate auxiliary variables
y_sim = z_sim+alpha*k_sim+(1-alpha)*l_sim;
i_sim = (Y*y_sim-C*c_sim-G*g_sim)/(delta*K);
p_sim = y_sim-l_sim;


