function model2=Simulate2(model2)
%%
T=model2.T;
sigma_u=model2.sigma_u;
sigma_epsilon=model2.sigma_epsilon;
sigma_xi=model2.sigma_xi;
rho=model2.rho;
Tr=model2.Tr;
Nsim=model2.Nsim;
kappa=model2.kappa;
rho_gamma=model2.rho_gamma;

ZTN=zeros(T,Nsim);

xt=ZTN;
it=ZTN;
gammat=ZTN;
MPS=ZTN;
gamma_rolling=ZTN;
%optimal forecast E(gamma_{t+1}|Y_t)
gammaplus_LL=ZTN;
%estimation uncertainty Var(gamma_{t+1}|Y_t)
sigmat_LL=ZTN;
%fed funds forecast error
ff_error=ZTN;

b1=zeros(101,Nsim);
b2=zeros(101,Nsim);

b1r=zeros(101,Nsim);
b2r=zeros(101,Nsim);

coeffs_ff=zeros(4,Nsim);
%%
for i=1:Nsim
%draw random shocks
ut=mvnrnd(0,1,T)*sigma_u;
epsilont=mvnrnd(0,1,T)*sigma_epsilon;
xit=mvnrnd(0,1,T)*sigma_xi;

%simulate underlying process
for t=2:T
    xt(t,i)=rho*xt(t-1,i)+epsilont(t);
    gammat(t,i)=rho_gamma*gammat(t-1,i)+xit(t);
    it(t,i)=gammat(t)*xt(t,i)+ut(t);  
end

%rolling regressions
for t=Tr+1:T
    gamma_rolling(t,i)=(xt(t-Tr:t,i)'*it(t-Tr:t,i))/(xt(t-Tr:t,i)'*xt(t-Tr:t,i));
end

%% rational maximum likelihood estimation

%start with gamma known to be zero
gammaplus_LL(1,i)=0;
sigmat_LL(1,i)=0;

%updating 
for t=2:T
    %MPS
    MPS(t,i)=it(t,i)-gammaplus_LL(t-1,i).*xt(t,i);
    %Kalman filter with overconfidence parameter kappa
    omegat=sigmat_LL(t-1,i)^2/(sigmat_LL(t-1,i)^2+sigma_u^2/xt(t,i)^2);
    gammaplus_LL(t,i)=gammaplus_LL(t-1,i)+MPS(t,i)*omegat/xt(t,i);
    sigmat_LL(t,i)=kappa*(sqrt(sigmat_LL(t-1,i)^2*(1-omegat)+sigma_xi^2));
end

%% fed funds forecast error regressions
% ff_error(2:end,i)=it(2:end,i)-gammaplus_LL(1:end-1,i).*xt(1:end-1,i)*rho;
% 
% stdx     = std(xt(1:end-1,i));
% Delta_gamma = gammaplus_LL(2:end-1,i)-gammaplus_LL(1:end-2,i);
% std_Delta_gamma = std(Delta_gamma);
% lhsff    = ff_error(3:end,i); 
% rhsff    = [ones(size(ff_error(3:end,i))) xt(2:end-1,i)/stdx Delta_gamma/std_Delta_gamma xt(2:end-1,i).*Delta_gamma/(stdx*std_Delta_gamma)]; 
% varff   = rhsff\lhsff;
% coeffs_ff(:,i)=varff;

ff_error(5:end,i)=it(5:end,i)-gammaplus_LL(1:end-4,i).*xt(1:end-4,i)*rho^4;

stdx     = std(xt(1:end-4,i));
Delta_gamma = gammaplus_LL(5:end-4,i)-gammaplus_LL(1:end-8,i);
std_Delta_gamma = std(Delta_gamma);
lhsff    = ff_error(9:end,i); 
%rhsff    = [ones(size(ff_error(9:end,i))) it(5:end-4,i) xt(5:end-4,i)/stdx Delta_gamma/std_Delta_gamma xt(5:end-4,i).*Delta_gamma/(stdx*std_Delta_gamma)]; 
rhsff    = [ones(size(ff_error(9:end,i))) xt(5:end-4,i)/stdx Delta_gamma/std_Delta_gamma xt(5:end-4,i).*Delta_gamma/(stdx*std_Delta_gamma)]; 
varff   = rhsff\lhsff;
coeffs_ff(:,i)=varff;
%%
%state-contingent impulse responses
weak=(xt(:,i)<0);
for k=1:101
h=k-1;
lhs    = gammaplus_LL(2+h:end,i); 
rhs    = [ones(size(gammaplus_LL(2:end-h,i))) MPS(2:end-h,i).*(1-weak(2:end-h)) MPS(2:end-h,i).*weak(2:end-h) weak(2:end-h) gammaplus_LL(1:end-h-1,i)]; 

varb   = rhs\lhs;
b1(k,i)=varb(2);
b2(k,i)=varb(3);

lhs_rolling    = gammaplus_LL(2+h:end,i); 
rhs_rolling    = [ones(size(gammaplus_LL(2:end-h,i))) MPS(2:end-h,i).*(1-weak(2:end-h)) MPS(2:end-h,i).*weak(2:end-h) weak(2:end-h) gamma_rolling(1:end-h-1,i)]; 
varbr   = rhs_rolling\lhs_rolling;
b1r(k,i)=varbr(2);
b2r(k,i)=varbr(3);
end



end
b1_mean=mean(b1')';
b1_u=prctile(b1',90)';
b1_d=prctile(b1',10)';

b2_mean=mean(b2')';
b2_u=prctile(b2',90)';
b2_d=prctile(b2',10)';

b1r_mean=mean(b1r')';
b1r_u=prctile(b1r',90)';
b1r_d=prctile(b1r',10)';

b2r_mean=mean(b2r')';
b2r_u=prctile(b2r',90)';
b2r_d=prctile(b2r',10)';

coeffs_ff_mean=mean(coeffs_ff')';

model2.xt=xt;
model2.it=it;
model2.gammat=gammat;
model2.gamma_rolling=gamma_rolling;
model2.gammaplus_LL=gammaplus_LL;
model2.MPS=MPS;
model2.b1=b1;
model2.b2=b2;
model2.b1_mean=b1_mean;
model2.b1_u=b1_u;
model2.b1_d=b1_d;
model2.b2_mean=b2_mean;
model2.b2_u=b2_u;
model2.b2_d=b2_d;

model2.b1r=b1r;
model2.b2r=b2r;
model2.b1r_mean=b1r_mean;
model2.b1r_u=b1r_u;
model2.b1r_d=b1r_d;
model2.b2r_mean=b2r_mean;
model2.b2r_u=b2r_u;
model2.b2r_d=b2r_d;

model2.coeffs_ff_mean=coeffs_ff_mean;
end