In my last post I set up an example dataset and then showed how use the Durbin-Watson statistic to test for serial correlation. Once we have recognized that serial correlation exists, we need to make the appropriate adjustments because OLS is no longer the best estimation technique.
We need to rely on GLS to perform our estimation. The key to this is to identify what the matrix Omega (Z) looks like. In the case where we have first order autocorrelation AR(1), the Omega (Z) matrix looks like:
[1 p p2 p3 ... pn-1]
|p 1 p p2 ... pn-2|
|p2 p 1 p ... pn-3|
|.. .. .. .. ... ... |
[pn-1 pn-2 pn-3 pn-4 ... 1 ]
where the p is first order auto-correlation coefficient -- meaning the correlation between sequential error terms is p (rho). I would use Greek letters but I cant.)
So, we can easily apply GLS, we just need to figure out what p (rho) is.
One way of estimating p (rho) is to use the Durban-Watson statistic. It turns out that the DW Statistic is approximately equal to:
2 * (1 - p). In our regression, we computed a DW statistic of 1.09, which would imply a p (rho) of 0.455.
But, there is another, more interesting method, called the Hildreth-Lu Procedure.
This approach uses a grid search to estimate the most 'likely' value for p (rho). In practice this might involve running 100 regressions, each using a different value of p (rho), such as 0.01, 0.02, 0.03, ...0.99. We would then select the regression that yielded the lowest sum of squared errors. Once we have narrowed down the range, we could run another 100 regressions to further refine our estimate. What a perfect task for F# where we can send off multiple regressions at once and make full use of our dual core or multi-core computers.
We first transform our original data to look like this:
Y* = X* + v
Where
Y* = Yt - pYt-1
Each entry in X*: Xkt = Xkt - pXkt-1
v = et - pet-1
We then run a regular OLS regression and calculate the sum of squared errors. What we need is a function that will take a value of p (rho), do the data transformation and then run the OLS regression. The function would then return of the sum of squared errors.
Here's the code:
type rssType = {rho : float; rss : float}
let HL_iteration (X0 : float [,] ) (y : float [,] ) (p : float) =
let n1 = n-1
let X2 = Array2D.zeroCreate<float> n1 k
let y2 = Array2D.zeroCreate<float> n1 1
for i = 0 to (n1-1) do
y2.[i, 0] <- y.[i+1, 0] - p * y.[i, 0]
X2.[i, 0] <- 1.0 - p
X2.[i, 1] <- X0.[i+1, 1] - p * X0.[i, 1]
X2.[i, 2] <- X0.[i+1, 2] - p * X0.[i, 2]
X2.[i, 3] <- X0.[i+1, 3] - p * X0.[i, 3]
X2.[i, 4] <- X0.[i+1, 4] - p * X0.[i, 4]
let ret, RSS = regress X2 y2
let retVal = {rho=p; rss=RSS}
retVal
I first create a type rssType that can hold both rho and the sum of squared errors. The function takes our X and y matrices as input, as well a value for p (rho). We then get the count of observations and reduce it by one because we lose the first observation in the process (the first observation has no lagged value. We then create two arrays, X2 and y2 that will hold our transformed data. We loop over the original data and do the transformation element-by-element. Once we have our transformed data, we pass it to our regular OLS regression method.
We return the values of p (rho) and the sum residual sum of squares (sum or squared errors). The idea here is that we can call this function multiple times for different p's (rho's) and get back the sum of squares. We can then pick the rho with the smallest sum of squares.
This is how we do that:
let HL_iterate (X : float [,] ) (y : float [,] ) (i0 : float) (i1 : float) =
let incr = (i1 - i0) / 10.0
let rssList = new List<rssType>()
let asyncList = new List<Async<rssType>>()
for i = 0 to 9 do
let (iret : Async<rssType>) = async {return HL_iteration X y ((i0 + ((float) i * incr)))}
asyncList.Add iret
let Result =
async { let! asyncWorkflow = Async.Parallel asyncList
return asyncWorkflow }
let rssList = Async.RunSynchronously Result
for aRSS : rssType in rssList do
printfn "RSS, rho: %10.5f %10.6E " aRSS.rss aRSS.rho
rssList
This function takes our original X and y matrices, and two end points i0 and i1. The plan here is to take the two end points for rho, say 0.01 and 0.99, and split them into 10 equally spaced values for rho. We then run a regression for each of the ten values and obtain the results. By narrowing then end points we can fine tune our grid search to obtain a more precise estimate of rho. The variable incr is just the range i1 - i0 divided by 10. I then create a list rssList of type rssType to hold my results.
Now, I want to send off all my ten regressions at once, and take advantage of F#'s handling of multi-threaded operations. So, I set up a list of function calls I want to run asynchronously. These will be held in my variable asyncList. I next create a loop with 10 iterations than creates the function call for each of my spaced out values of rho. The value of rho: (i0 + ((float) i * incr) is just the starting point plus i times the value of the i of the increment variable. I add each function call to my asyncList.
I then set up the workflow to process the function calls in parallel and only return when they are all complete (this is the magic of the let! statement).
let Result =
async { let! asyncWorkflow = Async.Parallel asyncList
return asyncWorkflow }
I return the values in the variable asyncWorkflow which is a list of rssType values. Note: this statement does not actually run the functions (regressions). You actually need the following statement:
let rssList = Async.RunSynchronously Result
to run the regressions and fill up the rssList of results.
I then create a main function to call my HL_iterate function:
let rssList = HL_iterate X y 0.01 0.99
for aRSS in rssList do
printfn "RSS, rho: %10.6E %10.5f " aRSS.rss aRSS.rho
let rssArr = Array.sortWith (fun (x : rssType) (y : rssType) -> if y.rss = x.rss then 0 elif y.rss < x.rss then 1 else -1 ) rssList
printfn "Lowest RSS, rho: %10.6E %10.5f " (rssArr.[0].rss) (rssArr.[0].rho)
---------------------------------------------------------------------------
The variable rssList has the results from the ten regressions. I then print them out. My next step is to find the regression with the lowest sum of squares. I do this by invoking the Array method sortWith. This method allows me to sort an array using an anonymous function that I provide.
(fun (x : rssType) (y : rssType) -> if y.rss = x.rss then 0 elif y.rss < x.rss then 1 else -1 )
My function takes the value of rss (one part of the rssType structure) and compares it to the next value. My comparison uses the '<' operator so that I can get a resulting array that is sorted in ascending order. Thus, it is easy to obtain the regression with the lowest sum of squares -- it ends up at the head of the list.
So, let's run it and see what we get:
RSS, rho: 5.018932E+009 0.01000
RSS, rho: 4.636517E+009 0.10800
RSS, rho: 4.348996E+009 0.20600
RSS, rho: 4.156232E+009 0.30400
RSS, rho: 4.058028E+009 0.40200
RSS, rho: 4.054117E+009 0.50000
RSS, rho: 4.144200E+009 0.59800
RSS, rho: 4.328156E+009 0.69600
RSS, rho: 4.606552E+009 0.79400
RSS, rho: 4.981314E+009 0.89200
Lowest RSS, rho: 4.054117E+009 0.50000
I just love to see the both processor cores spike! We see that the lowest rho is in the middle at 0.50. We can fine tune this by setting our end points to 0.5 +/- 0.1 and calling the HL_iterate function again:
HL_iterate X y (rssList.[0].rho - 0.1) (rssList.[0].rho+0.1)
let rssList2 = HL_iterate X y (rssArr.[0].rho - 0.1) (rssArr.[0].rho + 0.1)
for aRSS in rssList2 do
printfn "RSS, rho: %10.6E %10.5f " aRSS.rss aRSS.rho
let rssArr2 = Array.sortWith (fun (x : rssType) (y : rssType) -> if y.rss = x.rss then 0 elif y.rss < x.rss then 1 else -1 ) rssList2
printfn "Lowest RSS, rho: %10.6E %10.5f " (rssArr2.[0].rss) (rssArr2.[0].rho)
We now get our second round of results
RSS, rho: 4.059089E+009 0.40000
RSS, rho: 4.050250E+009 0.42000
RSS, rho: 4.045336E+009 0.44000
RSS, rho: 4.044344E+009 0.46000
RSS, rho: 4.047272E+009 0.48000
RSS, rho: 4.054117E+009 0.50000
RSS, rho: 4.064876E+009 0.52000
RSS, rho: 4.079547E+009 0.54000
RSS, rho: 4.098129E+009 0.56000
RSS, rho: 4.120619E+009 0.58000
Lowest RSS, rho: 4.044344E+009 0.46000
We now add a third iteration, where we adjust the endpoints to be 0.46 +/- 0.01:
let rssList3 = HL_iterate X y (rssArr2.[0].rho - 0.01) (rssArr2.[0].rho + 0.01)
for aRSS in rssList3 do
printfn "RSS, rho: %10.6E %10.5f " aRSS.rss aRSS.rho
let rssArr3 = Array.sortWith (fun (x : rssType) (y : rssType) -> if y.rss = x.rss then 0 elif y.rss < x.rss then 1 else -1 ) rssList3
printfn "Lowest RSS, rho: %10.6E %10.5f " (rssArr3.[0].rss) (rssArr3.[0].rho)
We now get our third round of results:
RSS, rho: 4.044350E+009 0.45000
RSS, rho: 4.044270E+009 0.45200
RSS, rho: 4.044230E+009 0.45400
RSS, rho: 4.044229E+009 0.45600
RSS, rho: 4.044267E+009 0.45800
RSS, rho: 4.044344E+009 0.46000
RSS, rho: 4.044460E+009 0.46200
RSS, rho: 4.044616E+009 0.46400
RSS, rho: 4.044811E+009 0.46600
RSS, rho: 4.045045E+009 0.46800
Lowest RSS, rho: 4.044229E+009 0.45600
So, I think I'll stop here. My best guess as rho is 0.456.
Now I'll just build my Z matrix (shown at the top) and apply my GLS estimation method.
This two-step procedure is commonly referred to as the Prais-Winsten procedure.
It turns out, given the matrix of the form:
[1 p p2 p3 ... pn-1]
|p 1 p p2 ... pn-2|
|p2 p 1 p ... pn-3|
|.. .. .. .. ... ... |
[pn-1 pn-2 pn-3 pn-4 ... 1 ]
that the inverse is easily computed. The inverse of Z is:
Z-1 = (1 / (1 - p2) * [ 1 -p 0 0 0 0]
|-p 1+p2 -p 0 0 0|
| 0 -p 1+p2 -p 0 0|
| .. .. .. .. .. ..|
| 0 0 0 -P 1+p2 -p|
[ 0 0 0 0 -p 1]
So for our program all we have to do is create Z-1 based on our estimate of p (rho) and then run our GLS regression
We now add the following to our program:
let rho = rssArr3.[0].rho
//----------------------------------------------------------------------
// Now Build ZI Matrix
//----------------------------------------------------------------------
let dd = (1.0 / (1.0 - Math.Pow(rho,2.0)))
let ZI = Array2D.init<float> n n (fun i j -> if i=j then dd * (1.0 + Math.Pow(rho,2.0)) elif (i=j+1 or j=i+1) then -1.0*rho*dd else 0.0)
ZI.[0,0] <- dd * 1.0
ZI.[n-1,n-1] <- dd * 1.0
let ret = regressGLS X y ZI
We first grab rho estimate from our last Hildreth-Liu iteration. Recall it should be 0.456. We then build a ZI matrix, which is Omega Inverse.
The variable dd is just the multiplier out in front of the matrix (= 1 / (1 - p2)).
ZI is an n x n matrix where the diagonal elements are 1 + p2 and the adjacent entries are merely -p. The diagonal corners are set to 1.0. All these elements are multiplied by dd.
I then call my regressGLS method, passing to it our data matrices and Omega inverse (ZI). The full code of RegressGLS is shown below:
--------------------------------------------------------------------------------
let regressGLS (x0 : float [,] ) (y : float [,] ) (ZI : float [,] ) =
let xp = transpose2 x0
//--------------------------------
// Now create X'ZI X
//--------------------------------
let xpx = matMult2 (matMult2 xp ZI) x0
let X = xpx
let XI, ret = matInverse X
let k = X.GetUpperBound 0
let ni = 1 + x0.GetUpperBound 0
for i = 0 to k do
//printfn "XI Matrix: %10.5f %10.5f " XI.[i,0] XI.[i,1]
printfn "XI Matrix: %15.10f %15.10f %15.10f %15.10f %15.10f " XI.[i,0] XI.[i,1] XI.[i,2] XI.[i,3] XI.[i,4]
//-----------------------------------------
//Now Test
//Multiply XI and X to calculate Identity Matrix
//-----------------------------------------
let XTest = Array2D.zeroCreate<float> (k+1) (k+1)
let XTest = matMult XI X XTest
for i = 0 to k do
//printfn "XTest Matrix: %10.5f %10.5f " XTest.[i,0] XTest.[i,1]
printfn "XTest Matrix: %10.5f %10.5f %10.5f %10.5f %10.5f " XTest.[i,0] XTest.[i,1] XTest.[i,2] XTest.[i,3] XTest.[i,4]
//let a0 = Array2D.zeroCreate<float> (k+1) ni
//let a1 = Array2D.zeroCreate<float> ni 1
let b = matMult2 (matMult2 (matMult2 XI xp) ZI) y
for i = 0 to k do
printfn "Beta Hat: %i %10.5f" i b.[i,0]
//---------------------------------------
// Calc Std Errors
//---------------------------------------
let yp = transpose2 y
let S1 = matMult2 (matMult2 yp ZI) y
let S2 = matMult2 (matMult2 (matMult2 (transpose2 b) xp) ZI) y
printfn "S1: %10.2f" S1.[0,0]
printfn "S2: %10.2f" S2.[0,0]
let sse = S1.[0,0] - S2.[0,0]
printfn "SSE: %10.2f" sse
let sigma2 = sse / float (ni - (k+1))
printfn "SE: %10.2f" sigma2
//--------------------------------------------------
//Calc Variance/Covariance Matrix of B Hat
//--------------------------------------------------
//---------------------------------------
// Calc T-Stats
//---------------------------------------
let bVar = Array2D.zeroCreate<float> (k+1) 1
let bSE = Array2D.zeroCreate<float> (k+1) 1
let TStat = Array2D.zeroCreate<float> (k+1) 1
for i = 0 to k do
bVar.[i,0] <- sigma2 * XI.[i,i]
bSE.[i,0] <- Math.Sqrt bVar.[i,0]
TStat.[i,0] <- b.[i,0] / bSE.[i,0]
printfn "B, Var, SE, TStat: %10.5f %10.5f %10.5f %10.5f" b.[i,0] bVar.[i,0] bSE.[i,0] TStat.[i,0]
--------------------------------------------------------------------------------
When we run the GLS regression here is what we get:
---------------------------------------------------------------------
Beta Hat: 0 31086.62449
Beta Hat: 1 -65.68368
Beta Hat: 2 739.92914
Beta Hat: 3 -85.11546
Beta Hat: 4 11.52357
S1: 932945003919.98
S2: 927832435762.29
SSE: 5112568157.69
SE: 11862107.09
B, Var, SE, TStat: 31086.62449 7459397.27394 2731.18972 11.38208
B, Var, SE, TStat: -65.68368 166.51318 12.90400 -5.09018
B, Var, SE, TStat: 739.92914 268.59157 16.38876 45.14856
B, Var, SE, TStat: -85.11546 414.70682 20.36435 -4.17963
B, Var, SE, TStat: 11.52357 20.11106 4.48454 2.56962
---------------------------------------------------------------------
My original straight OLS estimates are shown here:
---------------------------------------------------------------------
B, Var, SE, TStat: 29843.97538 3520606.31417 1876.32788 15.90552
B, Var, SE, TStat: -61.85082 74.88767 8.65377 -7.14727
B, Var, SE, TStat: 732.56913 126.54212 11.24909 65.12250
B, Var, SE, TStat: -68.85197 198.90621 14.10341 -4.88194
B, Var, SE, TStat: 10.10098 10.35856 3.21847 3.13844
---------------------------------------------------------------------
It turns out that the parameter estimates (Beta Hats) are pretty much the same between the two regressions. I suppose this is to be expected since both yield unbiased estimates of the parameters. The standard errors are higher in the GLS regression, but that is OK. The standard errors calculated under plain vanilla OLS are biased so they aren't any good anyway. The standard errors under GLS are unbiased and are the appropriate ones to use for hypothesis testing.
So, there you have F# & GLS!
PS. If anyone finds a problem, please let me know!