% Computes task 31 'Robustifying Muth Filter'
% Reference:  Hansen-Sargent "Robust Filtering and Control ...'
% Plots variance of the worst case shock as a function of parameters 
%  	and given filter. It is expected that the minimum occurs at 
%	the robust filter 
 
% Uses doublex8.m to compute optimal robust filter 
%      ssvar2.m    to compute steady state variance of the worst case shock 
% Input: none 
% Uses: task31reg.m (Optimal linear regulator specialized to the task) 
%	ssvar2.m    (Steady state Variance formula for worst shock) 
%       h.m         (forecast error generator vector function 
%	  	    h(alpha,theta,N1,N2,K,omega )
%	  omega can be a vector so the result is length(omega)x2 matrix

clear
% set parameters 
format long;

ci = i; % ci is now stores sqrt(-1); 

A = 1; 
C = 1; 

alpmat = linspace(1.2,1.78,30);
themat = [10.^(8+[0:-0.25:-7])]; % the theta values used in plots

nmat = length(themat);

		E = [0 1]; 
     Q=1; 
      
for m = 1:length(alpmat)      

for n = 1:nmat
   
     alpha = alpmat(m);
    
    	D = [alpha 0]; 

		theta = themat(n);

		sigma =-1/theta;
 
		% computing optimal filter with preference for robustness sigma 
      
     khat(m,n) = doublex8(A,C,D,E,Q,sigma); 

end

end

figure(1)

plot(log(themat),khat(length(alpmat),:))
%title('K(\theta)')
xlabel('log(\theta)')
ylabel('K(\theta)')
hold on
plot(log(themat),khat(length(alpmat),:),'*')

print -deps rmuth8.eps


figure(2)

surf(log(themat),alpmat,khat)
%title('K(\theta)')
view([2 -3 1.8])
xlabel('log(\theta)')
ylabel('\alpha')
zlabel('K(\theta)')

print -deps rmuth9.eps
