%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MATLAB code: Exercise 6.5
% File: exercise_6-5.zip
%
% BL estimation with estimators beta-widehat and beta-tilde.
% 
% BL model: Y_{t} = beta * eps_{t-1} * Y_{t-2} + eps_{t}, 
%           with eps_{t} ~ N(0,1) distributed.
%
% OUTPUT:
% CovCLS    = Coverage probability of beta-widehat estimator
% CovT      = Coverage probability of beta-tilde estimator
% LengthCLS = Average length 95% CI: beta-widehat estimator
% LengthT   = Average length 95% CI: beta-tilde estimator
%
% Reference:
% Giordano, F. and Vitale, C. (2003).
%   CLS asymptotic variance for a particular relevant bilinear time 
%   series model.
%   Statistical Methods & Applications, 12(2), 169-185.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n    = 500;
m    = 100;
beta = 0.5;
iter = 1000;
sum1 = 0;
sum2 = 0;
VY   = 1/(1-beta*beta);     % Var[Y_{t}]

for it=1:iter
   eps  = normrnd(0,1,n+m,1);
   y    = zeros(n+m,1);
   for i=3:n+m
       y(i) = beta*eps(i-1)*y(i-2)+eps(i);
   end

   data=y(m+1:m+n);

%   Start: CLS estimation    
    y2    = data.*data;
    y2lag = y2(3:n);
    yy    = y2(1:n-2);
    X     = [ones(size(y2lag)) yy];
    a     = X\y2lag; % a(1) is estimated LS variance of eps_{t} 

    num12  = 0;
    den    = 0;
    num113 = 0;
    num014 = 0;
    for i=1:n-4
       num12 = num12 + data(i+2)*data(i+1)*data(i); 
       den   = den + data(i)*data(i);                 
    end
    c0 = den/(n-4);
      
    bestgr = (num12/(a(1)*den));  
    cond1  = bestgr^8;       
    if (cond1 < (1/105))
      sum1             = sum1+1;          
      bgr(sum1,1)      = bestgr;
      biasgr(sum1,1)   = bestgr-beta;         
      varestgr(sum1,1) = (1+22*bestgr^2+9*bestgr^2-6*bestgr^2*VY)/(VY);
    end
       
    besttil = num12/(n*VY);       
    cond2   = besttil^8;
    if (cond2 < (1/105))
      sum2             = sum2+1;
      btil(sum2,1)     = besttil;
      biastil(sum2,1)  = besttil-beta;         
      termtil          = (1-3*besttil^4)*(1-15*besttil^6);
      vartilde(sum2,1) = (1-besttil^2)*(183*besttil^6+...
                         42*besttil^4+14*besttil^2+1)/termtil; 
    end   
end     
 
meanbgr      = mean(bgr);                    
stdbestgr    = mean(sqrt(varestgr))/sqrt(n);  

meanbtil   = mean(btil);                    
stdbesttil = mean(sqrt(vartilde))/sqrt(n);  

%  Percentage of times beta is in the 95% CI
   sum3 = 0;
   for i=1:sum1   
       lb1(i) = bgr(i,1) - 1.96* sqrt(varestgr(i,1))/sqrt(n);
       ub1(i) = bgr(i,1) + 1.96* sqrt(varestgr(i,1))/sqrt(n);
       lengthgr(i) = ub1(i)-lb1(i);
       if (beta>=lb1(i) && beta<=ub1(i)) 
           sum3 = sum3+1;
       end
   end
   CovCLS    = (sum3/sum1)*100;  
   LengthCLS = mean(lengthgr);   
   
   sum4 = 0;
   for i=1:sum2                     
       lb2(i) = btil(i,1) - 1.96* sqrt(vartilde(i,1))/sqrt(n);
       ub2(i) = btil(i,1) + 1.96* sqrt(vartilde(i,1))/sqrt(n);
       lengthtil(i)=ub2(i)-lb2(i);
       if (beta>=lb2(i) && beta<=ub2(i)) 
           sum4 = sum4+1;
       end
   end
   CovT    = (sum4/sum2)*100;  
   LengthT = mean(lengthtil);  
