format bank
load svardata
plot(dm,rmpy*diag([100 1 1 1]))
legend('r','M','P','Y');
lags=7;
nv=size(rmpy,2);
[By, Bx,u,xxi]=rfvar3(rmpy([1:302,331:end],:),7,ones(492,1),302,3,3);
dm(303:323)'
disp('These are omitted dates correspond to the Volcker reserve targeting period.')
disp('Hit any key to continue.')
pause
format short g
S=u'*u;
T=size(rmpy,1)-lags*2-28+nv+1;% nv+1 is number of dummy obs from lambda and mu non-zero.
                              % 21 is the number of observations deleted
idmat=[1 1 0 0;1 1 1 1;0 0 1 1;0 0 0 1]% Zeros of idmat define the identifying restrictions on A0
idx=find(idmat);
Arf=chol(fliplr(flipud(S/T)));
Arf=inv(fliplr(flipud(Arf)));
x=Arf(idx);
% flips are needed because chol computes *upper* triangular A s.t A'*A=S/T, whereas we need
% the corresponding *lower* triangular matrix.
[fh,xh,gh,Hh,itct,fc,rc]=csminwel('svarlh',x,.001*eye(9),[],1e-6,100,S,idmat,T)
disp('Hit any key to continue.')
pause
[fh1,xh1,gh,Hh,itct,fc,rc]=csminwel('svarlh',xh,Hh,[],1e-8,100,S,idmat,T)
[xh xh1]
Ah=zeros(4);
Ah(idx)=xh1
sh=zeros(4);
sh(idx)=sqrt(diag(Hh))
disp('The two matrices above together suggest that the distinction')
disp('between the first two equations (money supply and demand) is')
disp('weakly identified.  Why?')
disp('Hit any key to continue.')
pause
format long g
[fh1, T*sum(log(diag(chol(S/T))))+T*nv/2 ]
format short g
disp('Use the line above to evaluate, by the Schwarz criterion,')
disp('the accuracy of the structural restrictions.')
disp('Hit any key to continue.')
pause
resp=impulsdt(cat(3,Ah,prodt(Ah,By,2,1)),48);
% prodt is a tensor product.  See its program comments for its general
% interpretation.  Here, it multiplies each of the 4x4 By(j,:,:) matrices
% on the left by Ah.
plot(squeeze(resp(1,:,:))')
figure
plot(squeeze(resp(2,:,:))')
figure
plot(squeeze(resp(3,:,:))')
figure
plot(squeeze(resp(4,:,:))');
figure
plot(squeeze(resp(:,1,:))');
disp('These plots show that the responses to a policy shock are reasonable,')
disp('except for a severe price puzzle.')
disp('Label the graphs, save them, and use them to explain the basis for these claims')
disp('Hit any key to continue.')
pause
x0=x([1 3 4 6 8:end])
idmat=[1 1 0 0;0 1 0 0;0 0 1 1;0 0 0 1]
idx=find(idmat);
[fhr,xhr,gh,Hhr,itct,fc,rc]=csminwel('svarlh',x0,.001*eye(6),[],1e-8,100,S,idmat,T)
[fhr,xhr,gh,Hhr,itct,fc,rc]=csminwel('svarlh',x0,.001*eye(6),[],1e-8,100,S,idmat,T)
format long g
[fhr fh1]
format short g
disp('What is the difference between the restrictions we can check')
disp('with these two numbers and those we checked with the last pair')
disp('Hit any key to continue.')
pause
Ah=zeros(4);
idmat
idx=find(idmat);
Ah(idx)=xhr
resp=impulsdt(cat(3,Ah,prodt(Ah,By,2,1)),48);
w=[sin([0:.1:1]*pi/2)' cos([0:.1:1]*pi/2)'];
rresp=zeros([6,size(resp)]);
w=[sin([0:.2:1]*pi/2)' cos([0:.2:1]*pi/2)'];
for j=1:6
   rresp0=prodt(resp,[[w(j,1) -w(j,2);w(j,2) w(j,1)] zeros(2);zeros(2) eye(2)],2,2);
   rresp(j,:,:,:)=permute(rresp0,[1 3 2]);
end
% rresp is a 6x4x4x48 matrix.  Each rresp(i,:,:,:) is the full set of impulse responses 
% under a particular rotation of the first two equations.  Figuring out exactly
% what this means is part of the exercise.
figure;
plot(squeeze(rresp(:,1,1,:))');figure;
plot(squeeze(rresp(:,2,1,:))');figure;
plot(squeeze(rresp(:,3,1,:))');figure;
plot(squeeze(rresp(:,4,1,:))')
grid;
disp('This set of graphs might suggest that use of prior information')
disp('about the nature of responses to monetary contraction, together with allowing')
disp('for simultaneity, can resolve the price puzzle.  Explain.')
disp('Hit any key to continue.')
pause
figure;
plot(squeeze(rresp(:,1,2,:))');figure;
plot(squeeze(rresp(:,2,2,:))');figure;
plot(squeeze(rresp(:,3,2,:))');figure;
plot(squeeze(rresp(:,4,2,:))')
disp('Some of these graphs, or others you might generate, might')
disp('help with your discussion.')
disp('Hit any key to continue.') 
pause
figure;
plot(squeeze(rresp(:,1,3,:))');figure;
plot(squeeze(rresp(:,1,4,:))');figure;
plot(squeeze(rresp(3,:,1,:))')