/****************************************************************************
--- IPOL Estimates Quarterly real GDP From Annual Totals

---INPUT: Y - n*1 dimension matrix of annual observations
          Xstar  - 4n*2 dimension matrix of related qarterly observations
                   Two related series are retails sale and industrial 
                        production
---OUTPUT: Ystar - 4n*1 dimensional matrix of quarterly observations

---external PROCEDURES: proc find_a(a)

This program is written by Michael Boldin
---Original Source:
   Chow, Gregory C. and Lin An-loh: "Best Linear Unbiased Interpolation,
   Distribution, and Extrapolation of Time Series by Related Series",
   The Review of Economics and Statistics, Volume 53, Issue 4 (Nov.1971)
   372-375
---Notes: We estimate this model with ar(1) residuals corresponding to
   case 2 in Chow-Lin.
******************************************************************************/

library nlsys;  /*using nonlinear equations*/

n=rows(Y);                  /*sample size of annual series*/
convcrit=0.001;

Bmat = (eye(n).*.ones(1,4));    /*create the distribution matrix B*/
X = Bmat*Xstar;

a_0=.7;
norm = 10;
do while norm >= convcrit;

ii=1;
vstar = zeros(4*n,4*n);
/*now use the iteractive procedure described by Chow and Lin following 
equation
23*/

do while ii<=4*n;
   jj = 1;
   do while jj<=4*n;                     /*knowing a, build up vstar*/
      vstar[ii,jj] = a_0^(abs(ii-jj));
      jj=jj+1;
   endo;
   ii = ii+1;
   @ii;@
  endo;
  vmat = Bmat*vstar*Bmat';                /*knowing vstar and B, calculate 
v*/
  ivm=Inv(vmat);
  uhat = (eye(n)-X*INV(X'*ivm*X)*X'*ivm)*Y;       /*Chow-Lin equation (16)*/

  /*autocorrelation coefficient*/
  q_1 = (TRIMR(uhat,1,0)'*TRIMR(uhat,0,1))/(uhat'*uhat);
  "q_1"; q_1; @wait;@

  /*solve euqtions Chow-Lin equation (23)*/
  {a_1,f_of_a_1,jacb,rc}=NLSYS(&find_a,a_0);
  norm = abs(a_1-a_0); "norm" norm; @wait@
  a_0 = a_1;

endo;

/*now we can calculate Chow-Lin equation (15)*/
b = INV(X'*ivm*X)*X'*ivm*Y;
Ystar = Xstar*b + vstar*Bmat'*ivm*uhat;   /*Ystar is quarterly GDP*/
print Ystar;

/*procedure to solve for the quarterly ar coefficient from the annual ar 
coefficietnt,
the equation is thus different from that given in Chow-Lin equation (23)*/

proc find_a(a);
local num, den, fcn;
num = a^7+2*a^6+3*a^5+4*a^4+3*a^3+2*a^2+a;
den = q_1*(4+6*a+2*a^3+4*a^2);
fcn = num - den;
retp(fcn);
endp