%main31c
% modifies Chao and Sergei's main31b to do graphs
% 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 = [1,1.78]; % the alpha values used in plots

themat = [1e+8,7]; % the theta values used in plots

nmat = length(themat);

   
   % while is invoked to compute only the values we want, namely, first 
   % row and first column of d(alpha(i),theta(j)), where d is the variance
   
for n = 1:nmat
      
  for m = 1:length(alpmat)

    alpha = alpmat(m);
    D = [alpha 0]; 
    E = [0 1]; 
    Q=1; 
      
		% computing optimal filter for correctly specified model 
      
    kstar = 0.5*(sqrt(alpha^4+4*alpha^2)-alpha^2); 
    kstarmat(m,n) = kstar; 
    theta = themat(n);
    sigma =-1/theta;
 
		% computing optimal filter with preference for robustness sigma 
      
    khat = doublex8(A,C,D,E,Q,sigma);
    khatmat(m,n) = khat;
      
	  K = linspace(0.1,1,99); 
	  K = [K khat kstar];  % khat is appended so that we can get the minimum value exactly
	  del = zeros(length(K),1); % variance computed usual way
	  delfft = zeros(length(K),1); % variance computed with fft method

		% setting parameters for fft calculation

		nn=8;
		T=2^nn;
		omega = 2*pi.*[0:T-1]'/T;
		spect = zeros(length(K),T); 

		for i=1:length(K) 
 
  			[F] = task31reg(alpha, theta, khat); 
  			N1 = F(1); 
  			N2 = F(2); 
 
  			% compute steady state variance of the worst case shocks 
       % through ssvar2.m function-file 
           
       del(i) = ssvar2(N1,N2,alpha,K(i));
       delstar(m,n) = ssvar2(N1,N2,alpha,kstarmat(m,n));
  
  			% calculation of variance by fft
  
  			g = h(alpha,theta,N1,N2,K(i),omega);
  			
  			spect(i,:) = (g(:,1).*conj(g(:,1))+g(:,2).*conj(g(:,2)))';
		end; 
 
 		% select the minimum value of the variance, since khat is added to the end of K vector
 		% khat will always be chosen to be d(m,n), and j always equal the length of K
       
     d(m,n) = min(del); 
 
 
 		K = K(1:length(K)-2);% -2 to get rid of appended kstar and khat
 		delpt(:,n) = del(1:length(del)-2);
 

		% compare choices for the filters

		figure(n)
		plot(K',log(delpt(:,n)),'b-')
	%	title(['variance against K, \theta = ',num2str(theta)])
		xlabel('K')
     ylabel('\Delta')
     axis([0.1 1 0 4]) 
      
     % beautification 
		if m ~= length(alpmat)
         t1 = text(K(33),log(del(33)),['\alpha= ',num2str(alpmat(m))]);
		else
         t1 = text(K(40),log(del(40)),['\alpha= ',num2str(alpmat(m))]);
                 
		end

   hold on
   
   plot(khatmat(m,n),log(d(m,n)),'x')
   plot(kstarmat(m,n),log(delstar(m,n)),'+')
	 	
		if m == length(alpmat)
   		hold off
		end

   j = size(spect,1)-1;
   
   if m == 1 & n == 1
      figure(3)
      xlabel(['\omega  (\alpha = ',num2str(alpmat(m)),', \theta = ',num2str(themat(n)),')'])

   elseif m == 1 & n == 2
      figure(4)
      xlabel(['\omega  (\alpha = ',num2str(alpmat(m)),', \theta = ',num2str(themat(n)),')'])

   elseif m == 2 & n == 1
      figure(5)
      xlabel(['\omega  (\alpha = ',num2str(alpmat(m)),', \theta = ',num2str(themat(n)),')'])

   elseif m == 2 & n == 2
      figure(6)
   end
   
      spe = spect(j,:)';  % spectral decompostion for khat
      if n == 1
         plot(omega,spe,'-')
      else
         plot(omega,spe,'k-')
      end
      
      %title('spectral decomposition')
      if m == 2 & n == 2
         xlabel(['\omega  (\alpha = ',num2str(alpmat(m)),', \theta = ',num2str(themat(n)),')'])
      end
      
		if m ~= length(alpmat) & n == 1 %figure 3
         t2 = text(omega(90),spe(90),['\alpha= ',num2str(alpmat(m))]);
  		elseif m == length(alpmat) & n == 1 %figure 5
  			t2 = text(omega(75),spe(75),['\alpha= ',num2str(alpmat(m))]);
		end

		hold on

		 if n == 1
          plot(omega,spect(j+1,:)','g-.')
          %legend('k=khat','k=kstar')
      elseif n == 2 & m == 1
         plot(omega,spect(j+1,:)','r-.')
         %legend('k=khat','k=kstar')
      elseif n == 2 & m == 2
         plot(omega,spect(j+1,:)','k-.')
         %legend('k=khat','k=kstar')
         
      end

		axis([0 pi, 0 inf]) % symmetric, only draw half

 
	end % for m loop

end % for n loop
 
 hold off % hold off for figure (3)
 figure(3)
 xlabel(['\omega  (\alpha = ',num2str(alpmat(1)),', \theta = ',num2str(themat(1)),')']);
figure(4)
 xlabel(['\omega  (\alpha = ',num2str(alpmat(1)),', \theta = ',num2str(themat(2)),')']);

figure(5)
 xlabel(['\omega  (\alpha = ',num2str(alpmat(2)),', \theta = ',num2str(themat(1)),')']);


 
 % print the figures
 
 figure(1)
 print -deps rmuth1.eps
 
 figure(2)
 print -deps rmuth2.eps
 
 figure(3)
 print -deps rmuth3.eps
 
 figure(4)
 print -deps rmuth4.eps
 
 figure(6)
 print -deps rmuth6.eps
 
 
 
 
 

