PhillipsHannan - Multivariate Hannan Efficient

[TomDoan]

This is an updated version of a procedure originally written by Rob Trevor to do Hannan efficient with multiple RHS lag structures; particularly for use in Phillips' estimator for the cointegrating matrix with arbitrary multivariate correlation structure.

Code: Select all

*
* @PhillipsHannan( options )  B START END
* #<supp. card> list of (P) endogenous variables
* #<supp. card> list of (Q) predetermined variables
*
* This procedure uses Hannan's Efficient Estimator to: estimate a set of linear
* equations with stationary errors; and test linear restrictions on the parameter
* matrix.
*
* Parameters:
*   B      = Initial estimate of coefficient matrix (PxQ)
*   START  = First observation
*   END    = Last observation
*
* Options:
*   PHILLIPS/[NOPHILLIPS] Whether doing Phillips' estimation of cointegrating matrix
*   BANDSPEC/[NOBANDSPEC] (Valid only if PHILLIPS) Whether to calculate estimates
*      using Band Spectral technique instead of Hannan
*   TEST/[NOTEST] Whether testing a set of linear restrictions H*BVEC=R where
*      where BVEC is a column vector which contains the stacked rows of the matrix B                                               *
*   MEANCORR/[NOMEANCORR] Whether to remove a constant term
*   SEASONAL/[NOSEASONAL] Whether to remove seasonal terms
*   TREND/[NOTREND] Whether to remove a linear time trend
*   ORDINATE = Number of ordinates [depends upon length, seasonal]
*   WINDOW = [FLAT]/TENT
*   WIDTH = Window width [depends upon ordinates]
*   ITER = Number of iterations to update initial estimate [1]
*   OLS/[NOOLS] Whether iteration 1 estimates are constrained to be the
*      the OLS estimates
*   H  = A SxPQ array of S linear restrictions (H*BVEC=R) on the
*      vector of estimated coefficients BVEC (used with TEST option)
*   R  = A Sx1 vector of linear restrictions (H*BVEC=R) on the
*      vector of estimated coefficients BVEC (used with TEST option)
*
*  Output:  (GLOBAL arrays)
*     BVEC   = Coefficient matrix stacked into a vector (P*Qx1)
*     VVAR   = Covariance matrix of the parameter estimates, corresponding
*              to the vectorised array BVEC
*
*  Version 1.0
*  Copyright (c) Robert G Trevor August 1990, 1992
*  Trevor and Associates ACN 065 479 549
*  Email: rtrevor@alumni.princeton.edu
*
*  Reference:  Phillips, P.C.B. (1991).  "Error Correction and Long Run
*              Equilibria in Continuous Time." Econometrica, Vol. 59, pp. 967-980.
*
*  Used in:    Hall, V.B. and R.G. Trevor (1992).  "Long Run Equilibrium
*              Estimation And Inference: A Non-Parametric Application."  In
*              P.C.B. Phillips (ed), Models, Methods and Applications of
*              Econometrics: Essays in Honour of A.R. Bergstrom. Basil
*              Blackwell (1992).
*
*  Revision Schedule:
*   07/2005 Updated by Tom Doan, Estima, to version 6. Minor change to syntax
*           (converted H and R to options)
*   04/2008 Renamed PhillipsHannan to match file name
*
proc PhillipsHannan  b start end
type integer start end
type rectangular *b
*
option switch   phillips 0
option switch   bandspec 0
option switch   test     0
option switch   meancorr 0
option switch   seasonal 0
option switch   trend    0
option switch   ols      0
option int      ordinate
option int      width
option int      iter     1
option choice   window   1 flat tent
option rect     h
option vect     r
*
local index serlisty serlistx serlistz
local integer pp p q s nrow ncol startm endm startl endl
local integer nobs counter m k
local real scalefac delta x2 signif
local lvector serlabsy serlabsx
local rectangular cyx qyx trb cinvq u v rhs e
local symmetric  cyy qyy cee qee w vcor veqn
local rect cxx qxx
local vector omega beqn
local integer dobandspec
local integer i j jj iloop slead nords
local vect[series] d
local rect[complex] fyy fxx fyx fee feeinv
local rect[complex] cxupdate
*
declare symmetric vvar
declare vector bvec
*
local series seasons time
local equation equation
*
if phillips==0.and.bandspec==1
   {
   dis "*** WARNING BandSpec switch ignored when Phillips switch is off ***"
   compute dobandspec=0
   }
else
   compute dobandspec=bandspec
*
* Read in series names and calc number of series
*
enter(varying) serlisty
enter(varying) serlistx
inquire(matrix=serlisty) p
inquire(matrix=serlistx) q
*
* Set number of LHS variables
*
if phillips==1
   {
   compute pp=p+q
   }
else
   {
   compute pp=p
   }
*
* Check dimensions of coefficient matrix
*
compute nrow=%rows(b),ncol=%cols(b)
if nrow<>p.or.ncol<>q
   {
   dis "B must have as many rows as endogenous variables and"
   dis "  as many columns as predetermined variables"
   return
   }
*
* Check dimensions of restriction matrices
*
if test==1
   {
   *
   * Check dimensions of R
   *
   compute s=%rows(r),ncol=%cols(r)
   if ncol>1
      {
      dis "R must be a vector"
      return
      }
   *
   * Check dimensions of H
   *
   compute nrow=%rows(h),ncol=%cols(h)
   if ncol<>p*q
      {
      dis "H must have as many columns as the number of parameters (P*Q)"
      return
      }
   if nrow<>s
      {
      dis "H and R must have the same number of rows"
      return
      }
   }
else
   {
   *
   * Set up test that all coefficients are zero
   *
   compute S=P*Q
   }
*
* Store series labels
*
dimension serlabsy(p) serlabsx(q)
ewise serlabsy(i)=%l(serlisty(i))
ewise serlabsx(i)=%l(serlistx(i))
*
* Calculate start and end for input variables
*
inquire(reglist) startm endm
# serlisty serlistx
*
* Adjust for lags
*
if phillips==1
   {
   compute startm=startm+1
   }
if .not.%defined(start).and..not.%defined(end)
   {
   compute startl=startm
   compute endl=endm
   }
else
   {
   if start<startm.or.end>endm
      {
      dis "ERROR:  Not all series are defined over range"  START END
      dis "        Common range is                      "  STARTM ENDM
      return
      }
   compute startl=start
   compute endl=end
   }
*
* Define length of complex series
*
if %defined(ordinate)
   compute nords=ordinate
else
   compute nords=%freqsize(endl-startl+1)
*
* Set up data series
*
dim d(pp+q)
*
* Put Y series into 1 to PP
* Put X series into PP+1 to PP+Q
*
*
* set up deterministic components
*
if seasonal==1
   {
   seasonal seasons
   compute slead = (1:1 - 0:1) - 2    ; * calculate seasonal frequency - 2
   enter(varying) serlistz
   # constant  seasons{-slead to 0}
   }
else if meancorr==1
   {
   enter(varying) serlistz
   # constant
   }
else
   {
   dimension serlistz(0)
   }
if trend==1
   {
   set time = t
   enter(varying) serlistz
   # serlistz time
   }
*
* Regress variables against deterministic components. Stuff
* residuals into slots in "d"
*
do i=1,p
   linreg(noprint) serlisty(i) startl endl d(i)
   # serlistz
   label d(i)
   # "#"+serlabsy(i)
end do i
if phillips==1
   {
   do i=1,q
      linreg(noprint) serlistx(i) (startl-1) endl d(p+i)
      # serlistz
      set d(pp+i) startl endl = d(p+i){1}
      label d(pp+i)
      # "#l"+serlabsx(i)
      difference d(p+i) startl endl
   end do i
   }
else
   {
   do i=1,q
      linreg(noprint) serlistx(i) startl endl d(pp+i)
      # serlistz
      label d(pp+i)
      # "#"+serlabsx(i)
   end do i
   }
*
* Set up complex series
* FFT Y  are in     1            to           PP
*     X             PP+1                      PP+Q
* Periodograms YY:  PP+Q+1                    PP+Q+(PP*PP)
*              YX:  PP+Q+.5*PP(PP+1)+1        PP+Q+.5*PP(PP+1)+(PP*Q)
*              XX:  PP+Q+.5*PP(PP+1)+(PP*Q)+1 PP+Q+.5*PP(PP+1)+(PP*Q)+.5*Q(Q+1)
*
compute counter = pp + q + pp*(pp+1)/2 + pp*q + q*(q+1)/2 + 1
frequency counter nords
*
* Convert time domain data to frequency domain
*
rtoc startl endl 1
# d
# 1 to (pp+q)
*
* Take Fast Fourier Transforms
*
do i=1,(pp+q)
   fft i
end do i
*
* Calculate cross periodograms
*
compute nobs=endl-startl+1
compute scalefac = 1./(2.*%pi*nobs)   ; * RATS doesn't divide fft by sqrt(2*pi*t)
compute counter=pp+q
do i=1,pp
   do j=1,i
      compute counter = counter + 1
      cmult(scale=scalefac) i j / counter
   end do j
end do i
do i=1,pp
   do j=(pp+1),(pp+q)
      compute counter = counter + 1
      cmult(scale=scalefac) i j / counter
   end do j
end do i
do i=(pp+1),(pp+q)
   do j=(pp+1),i
      compute counter = counter + 1
      cmult(scale=scalefac) i j / counter
   end do j
end do i
if counter<>pp+q+pp*(pp+1)/2+pp*q+q*(q+1)/2
   {
   dis "error:  count error in loop:  location 1"
   return
   }
*
* Smooth cross periodograms
* This shifts series in by (PP+Q)
* Spectrums of YY are   1                      .5*PP(PP+1)
*              YX       .5*PP(PP+1)+1          .5*PP(PP+1)+(PP*Q)
*              XX       .5*PP(PP+1)+(PP*Q)+1   .5*PP(PP+1)+(PP*Q)+.5*Q(Q+1)
*
* Set up mask for frequency zero to over-ride default
*
cset (counter+1) = 1.
*
* Smooth
*
do i=(pp+q+1),counter
   window(type=window,width=width,mask=(counter+1)) i / (i-(pp+q))
   cmult (i-(pp+q)) (counter+1) / (i-(pp+q))
end do i
*
* Set up complex arrays
*
compute m = fix(nords/2.)
if phillips==1
   {
   dimension e(pp,p)
   ewise e(i,j)=(i==j)
   }
*
* For each iteration calculate estimates
*
do iloop=1,iter
   dimension cyy(pp,pp) qyy(pp,pp) cyx(pp,q) qyx(pp,q) cxx(q,q) qxx(q,q)
   dimension cee(pp,pp) qee(pp,pp) rhs(p,q) w(p*q,p*q)
   dimension fyy(pp,pp) fyx(pp,q) fxx(q,q) fee(pp,pp)
   dimension u(pp,pp) v(pp,pp)

   compute rhs = %const(0.0)
   compute w   = %const(0.0)
   do k=0,m
      if (dobandspec==1.and.k<>0).and.(ols==0.or.iloop<>1)
         {
         break ;* Exit loop to get Band Spectral estimates
         }
      *
      * For each frequency 0<=w<=Pi
      *
      compute counter=0
      *
      * Set up Fyy = 2.0*(Cyy -iQyy)
      *
      do i=1,pp
         do j=1,i
            compute counter=counter+1
            compute fyy(j,i) = 2*%z(k+1,counter)
            compute fyy(i,j) = %conjg(fyy(j,i))
         end do j
      end do i
      *
      * Set up Fyx = 2.0*(Cyx -iQyx)
      *
      do i=1,pp
         do j=1,q
            compute counter=counter+1
            compute fyx(i,j) = 2.0*%conjg(%z(k+1,counter))
         end do j
      end do i
      *
      * Set up Fxx = 2.0*(Cxx -iQxx)
      *
      do i=1,q
         do j=1,i
            compute counter=counter+1
            compute fxx(j,i) = 2.0*%z(k+1,counter)
            compute fxx(i,j) = %conjg(fxx(j,i))
         end do j
      end do i
      if counter<>pp*(pp+1)/2+pp*q+q*(q+1)/2
         {
         dis "ERROR:  COUNT ERROR IN LOOP:  LOCATION 2 -" K
         return
         }
      *
      * Calculate Fee = .5*(Cee -iQee)
      *
      if (ols==0.or.iloop<>1)
         {
         if phillips==1
            {
            compute trb = tr(e*b)
            compute fee = fyy + e*b*fxx*trb - fyx*trb - e*b*%cxadj(fyx)
            }
         else
            {
            compute trb = tr(b)
            compute fee = fyy + b*fxx*trb - fyx*trb - b*%cxadj(fyx)
            }
         }
      else
         {
         *
         * To get OLS estimates ignore spectrum of residuals - assume flat
         *
         ewise fee(i,j)=%if(i==j,1./%pi,0.0)
         }
      compute feeinv=%cxinv(fee)
      if k==0.or.k==m
         {
         if dobandspec==0.or.(ols==1.and.iloop==1)
            {
            compute delta=0.5/float(m)
            }
         else
            {
            compute delta=1.0
            }
         }
      else
         {
         compute delta=1.0/float(m)
         }
      if phillips==1
         {
         *
         * Replace Hannan's U and V by E'U and E'V
         *
         compute feeinv=tr(e)*feeinv
         }
      compute cxupdate=feeinv*fyx
      ewise rhs(i,j)=rhs(i,j)+delta*%real(cxupdate(i,j))
      if phillips==1
         {
         *
         * Replace Hannan's U and V by E'UE and E'VE
         *
         compute feeinv = feeinv*e
         }
      ewise u(i,j)=%real(feeinv(i,j))
      ewise v(i,j)=%imag(feeinv(i,j))
      ewise cxx(i,j)=%real(fxx(i,j))
      ewise qxx(i,j)=%imag(fxx(i,j))
      compute w = w + delta*(%kroneker(u,cxx) - %kroneker(v,qxx))
   end do k
   release cyy qyy cyx qyx cxx qxx cee qee trb cinvq u v
   *
   * Set OMEGA to the stacked rows of RHS
   *
   dimension omega(p*q)
   compute counter=0
   do i=1,p
      do j=1,q
         compute counter=counter+1
         compute omega(counter) = rhs(i,j)
      end do j
   end do i
   if counter<>(p*q)
      {
      dis "ERROR:  COUNT ERROR IN LOOP:  LOCATION 3"
      return
      }
   release rhs
   compute vvar = inv(w)
   compute bvec = vvar*omega
   compute vvar = vvar/float(nobs)
   if dobandspec==1
      {
      *
      * Degrees of Freedom scale factor
      *
      compute vvar = (2.*m)*vvar
      }
   release omega w
   *
   * Update B matrix
   *
   compute counter=0
   do i=1,p
      do j=1,q
         compute counter=counter+1
         compute b(i,j) = bvec(counter)
      end do j
   end do i
   if counter<>(p*q)
      {
      dis "ERROR:  COUNT ERROR IN LOOP:  LOCATION 4"
      return
      }
   *
   * Calculate correlations across parameter estimates
   *
   dimension vcor(p*q,p*q)
   do i=1,p*q
      do j=1,i
         compute vcor(i,j) = vvar(i,j)/sqrt(vvar(i,i)*vvar(j,j))
      end do j
   end do i
   *
   if iloop==1
      {
      dis
      if phillips==1
         {
         dis "                     PHILLIPS COINTEGRATION ESTIMATION"
         }
      if dobandspec==0
         {
         dis "                         HANNAN EFFICIENT ESTIMATES"
         }
      else
         {
         dis "                           BAND SPECTRAL ESTIMATES"
         }
      dis
      dis "Period:" %datelabel(startl) %datelabel(endl) "=>" #### (ENDL-STARTL+1) "obs"
      dis
      dis "Endogenous variables are:"
      write(noskip) serlabsy
      dis
      dis "Predetermined variables are:"
      write(noskip) serlabsx
      dis
      }
   dis "Iteration:" ILOOP
   dis "Parameter estimates: "
   write(noskip) b
   dis
   dis "Correlation matrix of parameter vector"
   write(noskip) vcor
   dis
   release vcor
end do iloop
release e
*
* Calculate test statistic and significance level
*
if test==1
   {
   compute x2 = %qform(inv(h*(vvar*(tr(h)))),r-h*bvec)
   }
else
   {
   compute x2 = %qform(inv(vvar),bvec)
   }
cdf(noprint) chisquared x2 s
compute signif=%signif
if test==1
   {
   dis "WALD TEST OF RESTRICTIONS H*BVEC=R"
   dis
   dis "H"
   write(noskip) h
   dis
   dis "R"
   write(noskip) r
   dis
   }
else
   {
   dis "WALD TEST OF RESTRICTION THAT ALL PARAMETERS ARE ZERO"
   dis
   }
dis "------------------- "
dis "Wald Stat | P-Value "
dis "------------------- "
dis   ######.## X2  "|"  #.#### SIGNIF
dis "------------------- "
dis
*
* Write out estimates equation by equation
*
dimension beqn(q) veqn(q,q)
*
* Change SERLISTY to transformed variables
*
enter(varying) serlisty
# d(1)
do i=2,p
   enter(varying) serlisty
   # serlisty d(i)
end do i
*
* Change SERLISTX to transformed variables
*
enter(varying) serlistx
# d(pp+1)
do j=2,q
   enter(varying) serlistx
   # serlistx d(pp+j)
end do j
equation(noconstant) equation serlisty(1)
# serlistx
*
* For each equation set up param vector and cov matrix
*
do i=1,p
   do j=1,q
      compute beqn(j) = bvec((q*(i-1))+j)
      do jj=1,j
         compute veqn(j,jj) = vvar((q*(i-1))+j,(q*(i-1))+jj)
      end do jj
   end do j
   dis
   dis
   if phillips==1
      {
      dis "                     PHILLIPS COINTEGRATION ESTIMATION"
      }
   if (ols==1.and.iter==1)
      {
      *
      * Turn on residual variance calculation
      *
      dis "                               OLS ESTIMATES"
      linreg(create,coeff=beqn,covmat=veqn,scale,equ=equation) $
               serlisty(i)   startl endl
      }
   else
      {
      if dobandspec==0
         {
         dis "                         HANNAN EFFICIENT ESTIMATES"
         }
      else
         {
         dis "                           BAND SPECTRAL ESTIMATES"
         }
      linreg(create,coeff=beqn,covmat=veqn,noscale,equ=equation) $
               serlisty(i)   startl endl
      }
end do i
*
release serlisty serlistx serlistz serlabsy serlabsx
release beqn veqn
end