' example subroutine for Fully-Modified OLS estimation (v1.0: 2009-04-28 - gts)
'
' *** note: this routine has undergone some testing and matches a variety
' of known results but is strictly use-at-your-own-risk ***
'
subroutine local fmols(series y, group g, scalar bw,  vector b, vector ses, sym cov)
'
' Arguments:
' 	y - dependent series
' 	g - cointegrating regressors
' 	bw - bandwidth (lag trunctation + 1)
' 	b - results coefficient vector
' 	ses - results standard errors
' 	cov - results covariance
'
' 1. only the Bartlett kernel is supported
'
' 2. estimator assumes a constant in the cointegrating equation and no additional 
' regressors in the cointegrating regressors equations
'
' 3. no df correction in covariance estimates
'
' example of code to call the routine on data in an EViews workfile:
' --------------------------------------------------------------------------------------
' smpl @all
' scalar bw = 7
' group g  x1 x2 x3
' vector b
' vector tstats
' sym cov
' call fmols(y, g, bw, b, ses, cov)
 ' vector tstats = @ediv(b, ses)
' --------------------------------------------------------------------------------------

' backup original sample
sample orig

' long run equation
!has_c  = 1
equation eq01.ls y c g
eq01.makeresid res1

' bartlett weights
vector(bw) bart
for !i=1 to bw
	bart(!i) = 1.0 - (!i / bw)
next

' form residuals
group gres res1 ' system residuals r
group gxres ' regressors residuals 
for !i=1 to g.count
	gres.add d(g(!i))
	gxres.add d(g(!i))
next

' get the sample information on the residual subset
smpl orig if @rnas(gres)=0
!obs = @obs(res1)

' estimate long-run covariance of the residuals using Bartlett kernel
!nvars = gres.count
do gres.cov(out=a) usscp
sym test1 = asscp
matrix test2 = test1
for !i=1 to bw
	group gresl gres ' form augmented block of resids: [ r  r(-!i) ]
	for !j=1 to gres.count
		gresl.add @lag(gres(!j), !i)
	next
	do gresl.cov(out=a) usscp
	sym c1 = asscp
	matrix rc1 = @subextract(c1, 1, !nvars+1, !nvars, 2*!nvars) ' extract upper block
	test1 = test1 + bart(!i)*(rc1 + @transpose(rc1))
	test2 = test2 + bart(!i)*rc1
next
test1 = test1 / !obs
test2 = test2 / !obs

' estimate long-run covariance of the residuals
sym omega = test1

' estimate the one-sided long-run covariance of the residuals
matrix gamma = test2

' get the swept matrix
sym iomega22 = @inverse(@subextract(omega, 2, 2))
matrix omega12 = @subextract(omega, 1, 2, 1)
matrix iomega12 = omega12*iomega22
matrix iomegaadj = iomega12*@transpose(omega12)

' get the adjustment factor
series uadj
group uadjg uadj
stomna(gxres, xres)
matrix temp = xres*@transpose(iomega12)
mtos(temp, uadjg)

' get bias adjustment gamma2+ = gamma2 - omega12 * Omega_22^{-1}*gamma2 (Phillips and Moon 5.17)
matrix gamma22 = @subextract(gamma, 2, 2)
matrix gamma2 = @subextract(gamma, 2, 1, @rows(gamma), 1) 
matrix adjust1 =  gamma2 - gamma22*@transpose(iomega12)
matrix lambdaplus1 = adjust1*!obs

' get the adjusted y+ = y - w12 Omega22^{-1} u2  (Hansen 1992 (13) and (14) p. 96) 
series yplus = y - uadj
series u1plus = res1 - uadj

' form the new group and compute the long-run covariance (one-sided)
' get the lambda+ (Hansen 1992, p. 95) - one sided long-run covariances of [ u1+, u2 ] using the corrected residual u1+ which
' avoids doing the larger matrix and then adjusting
group u u1plus gxres

' estimate long-run covariance of the residuals using Bartlett kernel
do u.cov(out=a) usscp
matrix test3 = asscp
for !i=1 to bw
	group gu u u1plus(-!i) 
	for !j=1 to gxres.count  ' form augmented block of resids: [ u  u(-!i) ]
		gu.add @lag(gxres(!j), !i)
	next
	do gu.cov(out=a) usscp
	sym c2 = asscp
	matrix rc2 = @subextract(c2, !nvars+1, 1, 2*!nvars, !nvars)
	test3 = test3 + bart(!i)*rc2
next
test3 = test3 / !obs

matrix lambda = test3
matrix lambdaplus = @subextract(lambda, 2, 1, @rows(lambda), 1)
lambdaplus = lambdaplus*!obs

' form the updated moment matrix for [ y+,  x ]
group data yplus c g 
do data.cov(out=out) usscp

' get the numerator xy+ with bias adjustment
vector(!nvars) xyadjust = 0
!xplace = !has_c + 1
matplace(xyadjust, lambdaplus, !xplace)
vector xyplus = @subextract(outsscp, 2, 1, @rows(outsscp), 1) - xyadjust

' get the coefficient vector
sym imoment = @inverse(@subextract(outsscp, 2, 2))
b = imoment* xyplus

' get the covariance matrix
scalar omega11= omega(1, 1) - iomegaadj(1, 1)
cov = omega11*imoment

' get the standard errors
ses = @sqrt(@getmaindiagonal(cov))

' restore the sample
smpl orig

endsub