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