The Mad Economist

While figuring out how to calculate the R-Square for a regression, I came across an F# method that is very cool.

    //---------------------------------------
    // Calc R-Squared
    //---------------------------------------
    let N = Array2D.zeroCreate<float> ni ni
    let I = Array2D.init<float> ni ni (fun i j -> if i=j then 1.0 else 0.0)
    let ix = Array2D.init<float> ni ni (fun i j -> 1.0/float ni)

    for i = 0 to (ni-1) do
        for j = 0 to (ni-1) do
            N.[i,j] <- I.[i,j] - ix.[i,j]   

    let sstV = matMult2 (matMult2 (transpose2 y) N) y

    let sst = sstV.[0,0]
    printfn "SST: %10.2f" sst
   
    let R2 = 1.0 - (sse / sst)
    printfn "R-Square: %10.6f" R2

    let RBar2 = 1.0 - (((float ni - 1.0) / (float ni - float k)) * (1.0 - R2))
    printfn "RBar-Square: %10.6f" RBar2

When instantiating a 2D array, you can use the method 'zeroCreate' that fills an n x k matrix with zeros.  But, if you want something other than zero you can use: 'init'.

Suppose I want to create an identity matrix that has ones along the main diagonal, but zeros everywhere else.  I can use this statement:

    let I = Array2D.init<float> ni ni (fun i j -> if i=j then 1.0 else 0.0)

to create it.  The input parameters are (1) the number of rows (2) the number of columns and (3) some anonymous function. 

The anonymous function is the tricky, but powerful part.  The syntax is: fun for function and then two integers -- I have used i & j but you can use anything as long as they are typed as int.  The function receives as input the row, column for each cell that is created.  Thus, I use the expression: if i=j then 1.0 else 0.0 to set ones along the main diagonal (meaning when i = j) and zeroes otherwise.  You can use any initializer function that uses the row/column address to set its value.

I use the statement:

    let ix = Array2D.init<float> ni ni (fun i j -> 1.0/float ni)

to create an ni x ni array where all cells are filled with the value (1/ni).  Note that I have to cast ni as a float before doing the division.

This is very cool.

So, why do I have to create these matrices?  Because, I need to calculate the total sum of squares.  The total sum of squares is:

TSS = y'Ny

where N is a n x n matrix that is filled with the value of (1 - (1/n)) on the diagonal and (1/n) otherwise.  What does this achieve?

Well, the total sum of squares equals the sum of (y - ybar)², meaning the sum of squared deviation of the y_i from the mean ybar.

So, y'Ny = y'Iy - y(1/n)y = Sum Y² - n(ybar)². which is what we want.

In the code above:

    let sstV = matMult2 (matMult2 (transpose2 y) N) y

sstV equals the product (matmult2) of the product of y transpose and N and y.  This line just does the two matrix multiplication operations in one line.

Then, R² is equal to: 1 - (SSE/SST).

The Rbar² (adjusted R²) equals 1 - ((n - 1)/(n - k)) * (1 - R²)

That's all there is to it.

So here are the results of the regression

Beta Hat: 0   36.86585
Beta Hat: 1    5.42270
Beta Hat: 2    9.14833
Beta Hat: 3    9.87322

Y, Hat, e: 0  290.68000  315.31145   24.63145
Y, Hat, e: 1  236.44000  255.74589   19.30589
Y, Hat, e: 2  397.78000  370.03833  -27.74167
Y, Hat, e: 3  209.13000  244.07443   34.94443
Y, Hat, e: 4  266.98000  273.71108    6.73108
Y, Hat, e: 5  359.83000  344.62716  -15.20284
Y, Hat, e: 6   56.08000   55.57818   -0.50182
Y, Hat, e: 7  299.26000  286.23658  -13.02342
Y, Hat, e: 8  254.52000  216.70835  -37.81165
Y, Hat, e: 9  496.51000  505.17856    8.66856

SSE:    4921.30
SE:     820.22

B, Var, SE, TStat:   36.86585  598.64049   24.46713    1.50675
B, Var, SE, TStat:    5.42270    0.64214    0.80134    6.76706
B, Var, SE, TStat:    9.14833    3.76546    1.94048    4.71447
B, Var, SE, TStat:    9.87322    3.30030    1.81667    5.43479

SST:  125033.78
R-Square:   0.960640
RBar-Square:   0.949395

----------------------------------------------------------------------------
----------------------------------------------------------------------------

Here's the complete code:

//-----------------------------------------
let regress =

    let X = xpx
    let XI, ret = matInverse X

    let k = X.GetUpperBound 0
   
    for i = 0 to k do
         printfn "XI Matrix: %10.5f %10.5f %10.5f %10.5f " XI.[i,0] XI.[i,1] XI.[i,2] XI.[i,3]

    //-----------------------------------------
    //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 %10.5f %10.5f " XTest.[i,0] XTest.[i,1] XTest.[i,2] XTest.[i,3]
   
    //let a0 = Array2D.zeroCreate<float> (k+1) ni
    //let a1 = Array2D.zeroCreate<float> ni 1
    let b = matMult2 (matMult2 XI xp) y

    for i = 0 to k do
         printfn "Beta Hat: %i %10.5f" i b.[i,0]
   
    //---------------------------------------
    //  Calc Std Errors
    //---------------------------------------

    // Calc Residual Vector = Y - Yhat
    // First Calc Y Hat = X BetaHat
   
    //let Y0 = Array2D.zeroCreate<float> ni 1
    let YHat = matMult2 x0 b
   
    let e = Array2D.zeroCreate<float> ni 1
    for i = 0 to (ni-1) do
        e.[i,0] <- YHat.[i,0] - y.[i,0]
       
    for i = 0 to (ni-1) do
         printfn "Y, Hat, e: %i %10.5f %10.5f %10.5f " i y.[i,0] YHat.[i,0] e.[i,0]

    //--------------------------------------------------
    //Calc SSE -> Sum of Squared Errors e'e
    //--------------------------------------------------
    //let zero2_e = Array2D.zeroCreate<float> 1 1
    let ssV = matMult2 (transpose2 e) e
   
    let sse = ssV.[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]
   

    //---------------------------------------
    // Calc R-Squared
    //---------------------------------------
    let N = Array2D.zeroCreate<float> ni ni
    let I = Array2D.init<float> ni ni (fun i j -> if i=j then 1.0 else 0.0)
    let ix = Array2D.init<float> ni ni (fun i j -> 1.0/float ni)

    for i = 0 to (ni-1) do
        for j = 0 to (ni-1) do
            N.[i,j] <- I.[i,j] - ix.[i,j]   

    let sstV = matMult2 (matMult2 (transpose2 y) N) y

    let sst = sstV.[0,0]
    printfn "SST: %10.2f" sst
   
    let R2 = 1.0 - (sse / sst)
    printfn "R-Square: %10.6f" R2

    let RBar2 = 1.0 - (((float ni - 1.0) / (float ni - float k)) * (1.0 - R2))
    printfn "RBar-Square: %10.6f" RBar2

    let userresp = Console.ReadLine()
    ()

Just to complete the OLS method, we need to add the F-Statistic:

    let FStat = (R2 / (1.0 - R2)) * ((float ni - float k) / (float k - 1.0))
    printfn "F Stat: %10.6f" FStat

That's all there is to it.