% SIMULS.M
%  Takes the output of SOLVES.M and simulates the state vector x(t) for
%  a recursive linear equilibrium model.  The stacked state vector
%  X(t) is formed, where X(t)'=[x(pt-p+1)', x(pt-p+2)', ... ,x(pt)'].
%  The law of motion for X(t) is
%
%    X(t+1) = Ah * X(t) + Ch * W(t+1),
%
%  where Ah=D^(-1)*F, and Ch=D^(-1)G, where
%          | I    0 0 ... 0 0  |       | 0  aop |
%          | -ao1 I 0 ... 0 0  |   F = |        |
%   D =    | 0 -ao2 0 ... 0 0  |       | 0   0  |
%          |                   |      | Cp  0  0 ...   0   |
%          |0 0 0 ... -aop-1 0 |      | 0  C1  0 ...   0   |
%                                 G = |                    |
%                                     |                    |
%                                     | 0   0  0 ...  Cp-1 |
%
%  and where W(t) is a vector white noise.  The output is stored in X.
%  Use SIMULC.M, SIMULI.M, SIMULK.M, SIMULG.M, SIMULB.M, and SIMULD.M
%  to generate paths for c, i, k, g, b, and d.
%  You will be asked to input the number of "years" that you want to simulate.
										
nT=input('State the number of "years" that you want to simulate,nT  ')
disp('The simulation will be of length nT X p')
n=max(size(ao1));
D=eye(n*p);
for j=1:p-1
nj=int2str(j);
eval(['D(j*n+1:(j+1)*n,(j-1)*n+1:j*n)=-ao',nj,';'])
end
F=zeros(n*p);
np=int2str(p);
eval(['F(1:n,n*(p-1)+1:n*p)=ao',np,';'])
[n1,n2]=size(C1);
G=zeros(n1*p,n2*p);
for j=2:p
nj=int2str(j);
eval(['G((j-1)*n1+1:j*n1,(j-1)*n2+1:j*n2)=C',nj,';'])
end
eval(['G(1:n1,1:n2)=C',np,';'])
Ah=D\F;
Ch=D\G;
[m1,m2]=size(Ch);
										
rand('normal')
u=rand(nT,m2);
CC=eye(n*p);
DD=zeros(n*p,m2);
xx0=zeros(1,n*p);
for j=1:p
nj=int2str(j);
eval(['xx0(1,n*(j-1)+1:n*j)=x0;'])
end
X=dlsim(Ah,Ch,CC,DD,u,xx0);
H=zeros(p,n*p)';
for j=1:p
nj=int2str(j);
eval(['H((j-1)*n+1:j*n,j)=sk',nj,''';'])
end
Y=X*H;YY=Y';
y=YY(:);
