The Mad Economist

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    p²    p³  ...  p^n-1]
|p    1    p     p²  ...  ^pn-2|
|p²   p    1     p   ...  p^n-3|
|..   ..   ..    ..  ...  ^... |
[p^n-1 p^n-2 p^n-3  p^n-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* = Y_t  - pY_t-1

Each entry in X*: X_kt = X_kt - pX_kt-1

v = e_t  - pe_t-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: