## August 2009 Blog Posts

Here is the chart for railroad traffic for the week ended 8/22/09
The data are for all non-intermodal rail carloads. I have been following this statistic as an indicator of the state of the economy. It appears that the trend started about 5 or 6 weeks ago of increasing rail traffic is holding to form. Total non-intermodal traffic stands at about 280,000 carloads versus 250,000 carloads (an increase of about 12 percent). While I was concerned a few weeks ago that the uptick was a blip, I am more confident now that we are seeing is the real deal. I should...

In my last post, I complained about the high cost of sugar quotas to US consumers, measuring about $1.8 billion a year. But, what does this have to do with a sugar shortage? In the last few months, the world sugar market has changed dramatically. Although historically the world price of sugar has been far below the US price, the gap has closed considerably. Below, I show the same price chart from my prior post:
Currently, the world price is 21.84 cents/pound and the US price is 26.45 cents/pound. The world price is actually above the target market price set by farm legislation...

I came across a story in the WSJ last week (WSJ August 13th) in which some major food companies were warning of a sugar shortage in the US. Excuse me. A sugar shortage? I'm not making this up. Here's the text of the letter:
August 5, 2009
The Honorable Thomas J. Vilsack
Secretary
U.S. Department of Agriculture
Jamie L. Whitten Federal Building
1400 Independence Avenue, S.W.
Washington, DC 20250
Dear Mr. Secretary:
The organizations and companies below urge you to increase the sugar import quota
immediately. Your experts forecast unprecedented shortages without prompt action.
According to USDA's World Agricultural Supply and Demand Estimates, the United
States will end the next fiscal year...

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|
|.. .. .. .. ... ...

So, In my quest to see where the economy stands, I came across some really enlightening data. Namely railroad traffic data published by the Association of American of Railroads. They report data on carloads, by type of product, on a weekly basis. I figure that looking at the weekly data can give a good idea of where the economy is in terms of the recession. So, here's the first chart:
This chart traces the weekly carload data from September 1, 2009 to the week ending August 8th 2009. The carload data are for all products, but exclude intermodal (trailers & containers) carloads. For reference,...

So, the next topic I want to take on is serial correlation.
For this example, I wanted to use some real data, so I pulled out my OLD copy of Pindyck and Rubinfeld (copyright 1976) and I borrowed one of their examples. They have an example that looks at coal consumption as a function of four variables: (1) iron & steel production, (2) electrical utility production, (3) the price of coal, and (4) the price of natural gas. The question here is to see how sensitive coal consumption is the relative price of coal versus natural gas. For example, it is...

In the general linear regression model, one of the assumptions is that the variance is constant across all observation points. Often times this assumption is not valid and variance is not uniform across observations. For example, if we were looking at how much people spend on housing based on their income, it probably is the case that the variance of the errors will be larger for high income folks than for poor folks. In this case, we call the errors "heteroskedastic" as opposed to being "homoskedastic". I just love trying to spell those words on an exam.
So, a more general...

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 =...

So we're finally ready to run a regression.
Let's first make up some data.
-----------------------------------------------------------
#light
open System
open System.Collections.Generic
open System.IO
open Xunit
open MathMod
open Microsoft.FSharp.Math
//---------------------------------------------------------
//Regression
//---------------------------------------------------------
let ni = 10
let k = 4
let x0 = Array2D.zeroCreate<float> ni k
let y = Array2D.zeroCreate<float> ni 1
y.[0,0] <- 290.68
y.[1,0] <- 236.44
y.[2,0] <- 397.78
y.[3,0] <- 209.13
y.[4,0] <- 266.98
y.[5,0] <- 359.83
y.[6,0] <- 56.08
y.[7,0] <- 299.26
y.[8,0] <- 254.52
y.[9,0] <- 496.51
x0.[0,0] <- 1.0
x0.[1,0] <- 1.0
x0.[2,0] <- 1.0
x0.[3,0] <- 1.0
x0.[4,0] <- 1.0
x0.[5,0] <- 1.0
x0.[6,0] <- 1.0
x0.[7,0] <- 1.0
x0.[8,0] <- 1.0
x0.[9,0] <- 1.0
x0.[0,1] <- 22.0
x0.[1,1] <- 4.0
x0.[2,1] <- 13.0
x0.[3,1] <- 3.0
x0.[4,1] <- 9.0
x0.[5,1] <- -13.0
x0.[6,1] <- -12.0
x0.[7,1] <- 7.0
x0.[8,1] <- 19.0
x0.[9,1] <- 21.0
x0.[0,2] <- 12.0
x0.[1,2] <- 14.0
x0.[2,2] <-...

In my last post, I showed how to take a matrix X and decompose it into two triangular matrices L (lower triangular) and U (upper triangular), such that X = LU
Can we get the inverse of X from L and U? Yes, indeedy.
The inverse of X equals the inverse of LU. One of the properties of matrices is that the inverse of the product of two matrices equals the product of the two inverses, where the ordering of the matrices is switched.
That is: if X-1 = (LU)-1 then X-1 = U-1L-1.
It turns out that because U and L are triangular, that...

Ok, so here is my F# code to perform a linear decomposition of a matrix:
First, let's make up some data:
open System
open System.Collections.Generic
open System.IO
open Xunit
open MathMod
open Microsoft.FSharp.Math
let X = Array2D.zeroCreate<float> 4 4
X.[0,0] <- 4.0
X.[1,0] <- 2.0
X.[2,0] <- 3.0
X.[3,0] <- 9.0
X.[0,1] <- 7.0
X.[1,1] <- 11.0
X.[2,1] <- 3.0
X.[3,1] <- -2.0
X.[0,2] <- -8.0
X.[1,2] <- 2.0
X.[2,2] <- -17.0
X.[3,2] <- 4.0
X.[0,3] <- 6.0
X.[1,3] <- 9.0
X.[2,3] <- 10.0
X.[3,3] <- 5.0
X is just a 4 x 4 array.
We next create a function that takes the X array as an input and does the decomposition, returning both L and U as the Lower and Upper matrices.
let deComp (X : float...