# Function the simulates the ACS with output growth rule and PLC solution
#
#--------------------------------------------------------------------------

function simulate_plc(params_in::Dict{String}{Float64},theta_in::Array{Float64,1},shocks_in::Dict{String}{Array{Float64,1}},init::Array{Float64,1},T::Int64,Tdrop::Int64)

#--------------------------------------------------------------------------
# Map Parameters
#--------------------------------------------------------------------------

tau    = params_in["tau"];
eta    = params_in["eta"];
psi1   = params_in["psi1"];
psi2   = params_in["psi2"];
rho_z  = params_in["rho_z"];
rho_g  = params_in["rho_g"];
rho_r  = params_in["rho_r"];
rho_d  = params_in["rho_d"];
chi_h  = params_in["chi_h"];
y_ss   = params_in["y_ss"];
R_ss   = params_in["R_ss"];
r      = params_in["r"];
pi_ss  = params_in["pi_ss"];
c_ss   = params_in["c_ss"];

#--------------------------------------------------------------------------
# Map Coefficients
#--------------------------------------------------------------------------

theta_pi_1::Array{Float64,1}, theta_yy_1::Array{Float64,1},theta_cc_1::Array{Float64,1}, theta_pi_2::Array{Float64,1}, theta_yy_2::Array{Float64,1},theta_cc_2::Array{Float64,1},d0::Float64,d1::Float64,d2::Float64,d3::Float64,d4::Float64,d5::Float64 = get_coefs_plc(theta_in,params_in);
#theta_pi_1, theta_yy_1,theta_cc_1, theta_pi_2, theta_yy_2,theta_cc_2,d0,d1,d2,d3,d4,d5 = get_coefs_plc(theta_in,params_in);

#--------------------------------------------------------------------------
#   Preallocate vectors
#--------------------------------------------------------------------------
TT = T + Tdrop;                                   # Total # of observations

ys  = zeros(TT,1);
cs  = zeros(TT,1);
hs  = zeros(TT,1);
ws  = zeros(TT,1);
pis = zeros(TT,1);
Rs  = zeros(TT,1);
zs  = zeros(TT,1);
gs  = zeros(TT,1);
ds  = zeros(TT,1);
ee1s= zeros(TT,1);
Rhats = zeros(TT,1);
yhats = zeros(TT,1);
pihats= zeros(TT,1);
chats= zeros(TT,1);

#--------------------------------------------------------------------------
# Map shocks
#--------------------------------------------------------------------------

ers = shocks_in["er"];
egs = shocks_in["eg"];
ezs = shocks_in["ez"];
eds = shocks_in["ed"];

#--------------------------------------------------------------------------
#   Main Simulation Code
#--------------------------------------------------------------------------

for tt = 1 : TT

    if tt == 1

        Rlag = init[1];
        ylag = init[2];
        g_lag = init[3];
        z_lag = init[4];
        d_lag = init[5];

    else

        Rlag = Rs[tt-1];
        ylag = ys[tt-1];
        g_lag = gs[tt-1];
        z_lag = zs[tt-1];
        d_lag = ds[tt-1];

    end

    Rhat_lag = log.(Rlag/R_ss);
    yhat_lag = log.(ylag/y_ss);

    gs[tt]   = rho_g*g_lag + egs[tt];
    zs[tt]   = rho_z*z_lag + ezs[tt];
    ds[tt]   = rho_d*d_lag + eds[tt];

    # basis function
    psi_use = [1, ers[tt],gs[tt],zs[tt],ds[tt],yhat_lag,Rhat_lag];

    Rbar = d0 + d1*ers[tt] + d2*gs[tt] + d3*zs[tt] + d4*ds[tt] + d5*yhat_lag;

    if Rhat_lag>Rbar

        # Approx. Policies at current point

        yhats[tt]  = dot(psi_use,theta_yy_1);
        pihats[tt] = dot(psi_use,theta_pi_1);

        chats[tt]  = yhats[tt]-gs[tt]

        # Transform to levels;
        pis[tt]   = pi_ss*exp(pihats[tt]);
        ys[tt]    = y_ss*exp(yhats[tt]);
        cs[tt]    = c_ss*exp(chats[tt]);

        Rhats[tt] = (1-rho_r)*psi1*pihats[tt] + (1-rho_r)*psi2*(yhats[tt]-yhat_lag + zs[tt]) + rho_r*Rhat_lag + ers[tt];

    else

        # Approx. Policies at current point
        yhats[tt]  = dot(psi_use,theta_yy_2);
        pihats[tt] = dot(psi_use,theta_pi_2);

        chats[tt]  = yhats[tt]-gs[tt]

        # Transform to levels;
        pis[tt]    = pi_ss*exp(pihats[tt]);
        ys[tt]     = y_ss*exp(yhats[tt]);
        cs[tt]    = c_ss*exp(chats[tt]);
        Rhats[tt] = log.(1/R_ss);

    end

    Rs[tt] = R_ss*exp(Rhats[tt]);

    # Additional objects
    hs[tt] = ys[tt];
    ws[tt] = chi_h*(ys[tt]^(1/eta))*cs[tt]^tau;

    # Check for errors


    if cs[tt]<0 || ys[tt]<0;
        println("Simulating negative consumption or output");
    end

    if isreal(pis[tt])==0;
        println("Simulating imaginary inflation")
    end

end

#--------------------------------------------------------------------------
# Trim simulation and output results
#--------------------------------------------------------------------------

outsim=Dict{String}{Array{Float64,1}}()
outsim["pi"]=pis[Tdrop+1:T+Tdrop]
outsim["R"]=Rs[Tdrop+1:T+Tdrop]
outsim["y"]=ys[Tdrop+1:T+Tdrop]
outsim["c"]=cs[Tdrop+1:T+Tdrop]
outsim["h"]=hs[Tdrop+1:T+Tdrop]
outsim["w"]=ws[Tdrop+1:T+Tdrop]
outsim["ee1"]=ee1s[Tdrop+1:T+Tdrop]
outsim["er"]=ers[Tdrop+1:T+Tdrop]
outsim["eg"]=egs[Tdrop+1:T+Tdrop]
outsim["ez"]=ezs[Tdrop+1:T+Tdrop]
outsim["ed"]=eds[Tdrop+1:T+Tdrop]
outsim["g"]=gs[Tdrop+1:T+Tdrop]
outsim["z"]=zs[Tdrop+1:T+Tdrop]
outsim["d"]=ds[Tdrop+1:T+Tdrop]
outsim["logy"]=log.(ys[Tdrop+1:T+Tdrop])
outsim["yhat"]=yhats[Tdrop+1:T+Tdrop]
outsim["pihat"]=pihats[Tdrop+1:T+Tdrop]
outsim["Rhat"]=Rhats[Tdrop+1:T+Tdrop]
outsim["Rlag"]=zeros(T)
outsim["ylag"]=zeros(T)

# Get lagged simulations of endogenous states %

if Tdrop==0
    outsim["Rlag"][1]= init[1];
    outsim["ylag"][1]= init[2];

    outsim["Rlag"][2:TT,1]=Rs[1:end-1];
    outsim["ylag"][2:TT,1]=ys[1:end-1];

else
    outsim["Rlag"]= Rs[Tdrop:end-1];
    outsim["ylag"]= ys[Tdrop:end-1];
end

return outsim

end
