%RBC_solve_sim - VERSION 2
%Solves and simulates RBC model with government spending
%Uses function RE_LIN_SOLVE.M
%Parameters eta and para are set in RBC_EVALUATE.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 = .25;
%para = 'nogov';
switch para
    case 'nogov'
        delta = .02;
        G = 0;
        gamma = 0;
        sigma_nu = 0;
    case 'gov'
        delta = .01;
        G = .9;
        gamma = .95;
        sigma_nu = .025;
end
        
%Parameters for simulation
T_init = 500;   %Number of initial values
T_sim = 284;    %Number of simulated observations  
N = 1000;       %Number of replications

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

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};

[PHI_U,PHI_S,B_S,PHI_v] = REE_lin_solve(Amat,Bmat,Cmat,PI);

%GENERATE SIMULATED DATA USING VECTORIZED METHOD
%Construct MAT for vectorized simulation
%Includes T_init initial values and 1 ending value
T = T_init+T_sim+1;
I = eye(3);
Ir = repmat(I,1,T);
Ic = mat2cell(Ir,size(I,1),repmat(size(I,2),1,T));
MAT_I = blkdiag(Ic{:});
PHIr = repmat(-PHI_S,1,T-1); 
PHIc = mat2cell(PHIr,size(PHI_S,1),repmat(size(PHI_S,2),1,T-1));
blkdiag(PHIc{:});
MAT_PHI = [zeros(3,3*T); blkdiag(PHIc{:}) zeros(3*(T-1),3)];
MAT = sparse(MAT_I+MAT_PHI);

%Run vectorized simulation
rng(0)
sig = [sigma_eps; sigma_nu];
innov_draw = randn(N*T,2).*repmat(sig',N*T,1);
innov_mat = reshape([zeros(T*N,1) innov_draw]',3*T,N);
sim_mat = MAT\innov_mat;
sim_repl = reshape(sim_mat,3,N*T)';
sim_data = zeros(T,N,3);
for j = 1:3
    sim_data(:,:,j) = reshape(sim_repl(:,j),T,N);
end
k_sim = sim_data(T_init+1:T-1,:,1);     %This is k_hat_t-1
z_sim = sim_data(T_init+1:T-1,:,2);
g_sim = sim_data(T_init+1:T-1,:,3);

lam_sim = k_sim*PHI_U(1)+z_sim*PHI_U(2)+g_sim*PHI_U(3);

c_sim = lam_sim*PHI_v(1,1)+k_sim*PHI_v(1,2)+z_sim*PHI_v(1,3)+g_sim*PHI_v(1,4);
l_sim = lam_sim*PHI_v(2,1)+k_sim*PHI_v(2,2)+z_sim*PHI_v(2,3)+g_sim*PHI_v(2,4);

%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);
w_sim = c_sim-l_sim/eta;
p_sim = y_sim-l_sim;
