% SIMPULS.M
%  Takes the output of SOLVES.M and creates two sorts of impulse response
%  functions for a periodic 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 and Ch are defined as in simul.m. and where W(t) is a vector
%  white noise.  The impulse response function of this system is obtained
%  and used to construct the two impulse response functions for the periodic
%  model.  The impulse response OF y(tp-p+k) (i.e., y in a particular
%  season) is put in d.  The impulse response TO w(pt-p+k) (i.e., to an
%  innovation in a particular season) is put in g.
%
%  The program prompts the user for the innovation number.  A few lines
%  at the end of the program are currently set up for p=4.
%  The user can edit the lines near the end of the program to
%  handle cases other than p=4.  The value of p must be in memory.
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)=si',nj,''';'])    % edit this line to change
end                                         % series
D=zeros(p,n*p);
HH=H';
%    ps is the number of shocks times p.
ps=nw*p;
for j=1:ps
nj=int2str(j);
eval(['Y',nj,'=dimpulse(Ah,Ch,HH,D,j,40);'])  % edit this to change number
end                                            % of terms evaluated
										
% nw is the number of shocks
i=input('give the shock number whose impulse response you want  ')
for j=1:p
nj=int2str(j);
ni=int2str(i);
eval(['d=Y',ni,'(:,j)'])
										
for k=nw+i:nw:p*n2-i+1
nk=int2str(k);
eval(['d=[Y',nk,'(:,j), d];'])
end
dd=d';
eval(['d',nj,'=dd(:);'])
end
%  d1 is the impulse response of first quarter y to current and past shocks,
%  d2 is the impulse response of second quarter y to current and past shocks,
%  and so on.  (Note: these are shifted by different amounts at the
%  beginning, and so must be adjusted)
										
%  Here's an  adjustment to align them that works when p=4:
%  You must edit this line if p is not 4.
%  We put the impulse response functions in d.
N=max(size(d1));
d=[d1(2*p-1:N), d2(2*p-2:N-1), d3(2*p-3:N-2), d4(2*p-4:N-3)];
Nd=max(size(d));
d=d(2:Nd,:);
%  Now compute impulse response for 'other' definition of the impulse
%  response, i.e., the response of y TO w(pt-p+k), the shock in a particular
%  period
										
for j=1:p
nj=int2str(j);
nh=int2str(i+(j-1)*nw);
eval(['Z=Y',nh,';'])
e=Z';
ee=e(:);
eval(['h',nj,'=ee;'])
end
										
%  Note:  These aren't aligned.  Here's an adjustment that aligns
%  them when p=4.
										
h=[h1(5:N-3),h2(6:N-2),h3(7:N-1),h4(8:N)];

plot(d),title('response of OF y(pt-p+k)'),pause
plot(h),title('response of y TO w(pt-p+k)'),pause
