
 %%  Estimatig a Sticky Price Model by Sbordone(2002) %%  
 %%  July 2003,  written by Ryo Kato  %%
 
 %% You need to load DATA prn file.
 %% Data file:  X = [ du dp  p-u ]
 %% X = 143(nearly 30 years) x 3 matrix, dULC, inflation, p/ULC
 %% Data source: BEA. See Hirata and Kato (2004) for details
 
 %% Load data 
 lag = 1;
    q   = 3 ;   
    load XR3.prn; 
    X = XR3 ;
    X = X(41:length(X),:) ;   %% Adjusting samples       
    MX = mean(X) ; 
    x1 = ones(length(X),1).*MX(1) ;
    x2 = ones(length(X),1).*MX(2) ;
    x3 = ones(length(X),1).*MX(3) ;
    XA = [x1  x2  x3 ];
    X = X -XA ;
    
 %% Preparing data for VAR(1) estimation
    Xt = X(2:length(X),:);
    X1 = X(1:length(X)-1,:);
    
 %% Creat constant terms   
   o1 = ones(length(X)-1,1);
  
   ZX =  [o1  X1];
   Y1 = Xt(:,1);
   Y2 = Xt(:,2);
   Y3 = Xt(:,3);
   
 %% Estimation of a VAR(1)  
   b1 = inv(ZX'*ZX)*ZX'*Y1 ;
   b2 = inv(ZX'*ZX)*ZX'*Y2 ; 
   b3 = inv(ZX'*ZX)*ZX'*Y3 ;
   
   bb = [b1  b2  b3 ];
   
   A = bb(2:4,:) ;
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
 Vs = []   ;
 for i = 1: 500;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
%% Compute various theoretical values for
%% p/ULC under various values of alpha0
  alpha0 = 27+i*0.002;   %% searching for optimal alpha0
  alpha1  = 1.00;
   b = 1 + alpha1 + 1/alpha0 ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
%% Get lamdas 
  lamda1 = (b - (b^2-4*alpha1)^0.5)/(2*alpha1);
  lamda2 = (b + (b^2-4*alpha1)^0.5)/(2*alpha1);
  lamda = [lamda1 ; lamda2  ] ;
   LD = (1/lamda2)*eye(q) ;
 %% Computing expected ULCs  
 Eigvalues = eig(LD*A) ;
 %%  Construtinc Eu
   J = inv(eye(q) - LD*A );
   EZ  = J*Xt' ;
   Eu  = EZ(1,:)';
   T3  = (1-lamda1)*Eu;
 %%  Arrange each terms
   T1  = lamda1*X1(:,3);
   T2  = Xt(:,1);
 PUstar = T1 - T2 + T3 ;  %% a theoretical value
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 PUreal = Xt(:,3) ;       %% actual data
 %% Calculate Variance of error
 %% Vars =100000 * var(PUreal-PUstar) ;
  d = PUreal-PUstar ;    %% residuals
  Vars = 100000 * mean((d-mean(d))'*(d-mean(d))) ; 
  Vs(i,1:2) = [Vars  alpha0];
  V1 = Vs(:,1);
  V2 = Vs(:,2);
  Min1 = min(V1) ;  %% minimized variance 
      
end
 %%  See where the minimun is!!
 
 C = [PUreal  PUstar];
 
   figure(1)
     subplot(2,1,1);
       plot(  V2,V1 );
       ylabel('variance')
       xlabel('alpha0');
     subplot(2, 1,2);
       plot( C );
       ylabel('Actual, Est with last alpha0')
       xlabel('time');
     
  
 
 
 
 
 
  