%This function simulates time series form the model
% output is in levels

% modified: KM 22/06/2005
% modified: CK 25/08/2005
function [outp,hrs,cons]=simulation(hor,inidrop);

[PAR_alpha, PAR_beta, PAR_delta, PAR_psi, PAR_gamma, PAR_g, PAR_mu,PAR_Lbar, PAR_taulbar, PAR_taux, PAR_rho, PAR_sigmatau, PAR_sigmaz]=parameters;
[C,K,L,Y,Cp,Kp,Lp,Yp] = SSfinder;

simhorizon=hor+inidrop;
load;
% x = [VAR_K VAR_X VAR_tau];
% y = [VAR_Y VAR_C VAR_L VAR_PROD];

     e=randn(simhorizon,2);
     e(:,1)=e(:,1)*PAR_sigmaz; 
     e(:,2)=e(:,2)*PAR_sigmatau/PAR_taulbar;
     
     eta=[0 0; 1 0; 0 1];
     x(:,1)=[0; 0; 0];
     z0=1;
     z(1)=z0;
     j=1;
     while j<simhorizon+1;
         xt=x(:,j);
         yt(:,j)=gx*xt;
         xold=xt;
         x(:,j+1)=hx*xold+eta*[e(j,1);e(j,2)];
         z(j+1)=z0*exp(PAR_mu+e(j,1));
         z0=z(j+1);
         j=j+1;
     end;
outp=z(1,inidrop:(hor+inidrop))'.*exp(Y+yt(1,inidrop:(hor+inidrop))');
cons=z(1,inidrop:(hor+inidrop))'.*exp(C+yt(2,inidrop:(hor+inidrop))');
hrs=exp(L+yt(3,inidrop:(hor+inidrop))');
%plot(z)
