function results = aggshocksim_fnc_RE(N,T,burn,Nsim,shocks,param,shock)

rho1 = param(2);
rho2 = param(3);
sigw = param(4);
sigv = param(5);

shockper=25+burn;

if shock==1
    shocks.w(1,shockper,:)=sigw*0.018;
end

for simct=1:Nsim

s=zeros(1,T+burn); 
stt=NaN(N,T); 
sttm1=NaN(N,T); 
stm1t=NaN(N,T); 
stm1tm1=NaN(N,T);

for t=3:T+burn
    s(1,t) = rho1*s(1,t-1) + rho2*s(1,t-2) + sigw*shocks.w(1,t,simct);
end

y = s + sigv*shocks.v;
AA = [rho1 rho2; 1 0]; 
BB = [sigw;0];
CC = [1 0];
DD = sigv;
s0=zeros(2,1); 
P0=eye(2);

% Run Kalman filter 
for n=1:N
    t=1;
        %predict
        sbm1 = AA*s0;
        Pbm1 = AA*P0*AA' + BB*BB';
        %update
        z  = y(n,t);
        K  = Pbm1*CC'/(CC*Pbm1*CC'+DD*DD');
        sb = sbm1 + K*(z-CC*sbm1);
        Pb = (eye(2)-K*CC)*Pbm1;
        stt(n,t)     = sb(1);
        stm1t(n,t)   = sb(2);
        sttm1(n,t)   = sbm1(1);
        stm1tm1(n,t) = sbm1(2);
        
    for t=2:T+burn
        %predict
         sbm1 = AA*sb;
         Pbm1 = AA*Pb*AA' + BB*BB';

        %update
        z  = y(n,t);
        K  = Pbm1*CC'/(CC*Pbm1*CC'+DD*DD');
        sb = sbm1 + K*(z-CC*sbm1);
        Pb = (eye(2)-K*CC)*Pbm1;
        stt(n,t)     = sb(1);
        stm1t(n,t)   = sb(2);
        sttm1(n,t)   = sbm1(1);
        stm1tm1(n,t) = sbm1(2);
    end 
end
% Remove burn-in 
ftt     = stt(:,burn+1:end);
ftm1t   = stm1t(:,burn+1:end);
fttm1   = sttm1(:,burn+1:end);
ftm1tm1 = stm1tm1(:,burn+1:end);
s       = s(burn+1:end);

%Compute one-step ahead consensus forecast
ftp1t = rho1*ftt   + rho2*ftm1t;
ftp2t = rho1*ftp1t + rho2*ftt;
ftp3t = rho1*ftp2t + rho2*ftp1t;
    
ftp1tm1 = rho1*fttm1   + rho2*ftm1tm1;
ftp2tm1 = rho1*ftp1tm1 + rho2*fttm1;
ftp3tm1 = rho1*ftp2tm1 + rho2*ftp1tm1;

ftt_mn   = mean(ftt,1);
ftp1t_mn = mean(ftp1t,1);
ftp2t_mn = mean(ftp2t,1);
ftp3t_mn = mean(ftp3t,1);

% Compute results
fe(:,simct) = s(1,2:end) - ftp1t_mn(1:end-1);
end

% Store results
results.fe=mean(fe,2);

end