## July 2009 Blog Posts

Well, it's now time to bite off the big enchilada -- matrix inversion.
So, what method to use. I first thought I would try the Gauss-Jordan method, but then when I picked up my copy of Numerical Recipes in C, I read the following:
"For inverting a matrix, Gauss-Jordan elimination is about as efficient as any other method. For solving sets of linear equations, Gauss-Jordan elimination produces both the solution of the equations for one or more right-hand side vectors b, and also the matrix inverse A. However, its principal weaknesses are (i) that is requires the right-hand sides to be stored...

Ok, so I'm poking around arrays in F#, and I come across a type called 'matrix.' What the devil is this?
You can create a matrix like this:
let aMatrix = matrix [ [ 1.0; 7.0 ];
[ 1.0; 3.0 ];
[ 1.0; 11.0 ];
[ 1.0; 9.0 ];
[ 1.0; 6.0 ] ]
Of course, there some bugs. If I run this in F# Interactive, I get:
"Tests.fs(10,15): error FS0039: The value or constructor 'matrix' is not defined. A construct with this name was found in FSharp.PowerPack.dll, which contains some modules and types that were implicitly referenced in some previous versions of F#....

Well, I have finally decided to write some econometrics routines in F#. I first did this type of stuff in FORTRAN when I was a grad student, so I figured I'd try to see if I could still figure it all out.
But, first things first.
We need routines to Transpose a matrix and to Multiply two matrices.
My transpose method is this:
--------------------------------------------------------------
//
// Let's make up some data
//
let n = 6
let k = 3
// Create data array - n x k
let x0 = Array2D.zeroCreate<float> n k
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.[0,1] <- 22.0
x0.[1,1] <- 4.0
x0.[2,1] <- 13.0
x0.[3,1]...

Well, against my better judgment, I decided to write a retirement calculator. The idea here is to calculate my monthly living allowance when I retire. How depressing. So here goes...
The problem is to take certain assumptions:
Current Age: let age = 45.0
Retirement Age: let retirementAge = 65.0
Current Wealth: let currentWealth = 200000.0
Rate of Return: let rr = 0.05
Age at Death: let death = 85.0
Annual Savings: let savings = 20000.0
Inflation: let pi = 0.02
Real RoR: let ir = rr - pi
And calculate the amount we can consume each month of retirement, leaving exactly zero behind. We have no heirs, except possibly the University of California. ...

Last post we wrote a heap sort routine.
Let's now take it for a road test:
[<Fact>]
let HeapSortTest() =
let testArray = [|11.5; 8.5; 9.5; 28.0; 25.5; 15.0; 12.0; 11.0; 9.0; 7.5|]
let xArray = MathMod.heapSort testArray testArray.Length
let yArray = Array.sort testArray
for i = 0 to testArray.Length-1 do
Assert.InRange(xArray.[i],(yArray.[i]-0.0001),(yArray.[i]+0.0001))
Again, I use xUnit to test whether the array is being sorted properly. I use the Array method '.sort' to build a comparison array, against which I can compare my results. Of course, I could just used the built-in sort method in F# for this task, but then I would never...

Today's addition to my math library is a heap sort function.
I often run into a problem with collections, where I have a collection of class objects that I need to access based on one of their characteristics. For example, suppose I have a collection of 10 automobile objects and I wan them sorted by their MPG.
My strategy here is this:
(1) Create an array of floating point numbers (the MPG for each car)
(2) Create a parallel array of integers that represent each car's index in the collection.
(3) Sort the floating point array by MPG, but use the auxiliary integer array to...

My question is: who cares?

Can anyone be so stupid? Curious.

I came across a cool tool yesterday for testing F# programs called xUnit and I thought I should give it a try, so here goes.
The first step: Get it and install it.
The latest version 1.5 is here: http://xunit.codeplex.com/Release/ProjectReleases.aspx?ReleaseId=28060
Don't worry about any documentation there is none -- you've just got to wing it. I downloaded the release and unzipped it:
You don't need to do any more that that.
The next thing I did was add a reference to the xunit.dll file:
Now for our first test. In my last post, I showed how to code a mortgage calculator. I now want to write a test...

Simple Problemo:
What's my monthly payment for a $200,000 fixed rate 30-year mortgage at 6.00%?
let mortgage = 200000.0
let n = 360
let i = 0.06
let ones = List.init n (fun one -> 1.0)
let pv0 = presentValue ones (i/12.0)
let levelPayment = mortgage / pv0
That's all there is to it. We first set up a mortgage of $200,000, the number of periods (months) of 360 and the interest rate of 6.00%. Next, we create a vector of $1.00 for n (360) periods.
The List function
List.init n (fun one -> 1.0)
creates a List with n members and a value equal to the result of some initializer...

Today's problem is:
Given an investment stream how to I calculate the internal rate of return (IRR)? For all of you home gamers out there, the internal rate of return is the discount rate that when applied to the investment stream yields a zero present value.
Take for example the following:
You have to pay $1,000 for an investment that returns the following cash flow over its life:
End of Year 1: $100.00
End of Year 2: $120.00
End of Year 3: $125.00
End of Year 4: $150.00
End of Year 5: $180.00
End of Year 6: $200.00
End of Year 7: $220.00
End of Year 8: $250.00
End of Year 9: ...

Today's problem is to calculate the present value of an investment stream.
Suppose you had to pay $1,000 for an investment that returns the following cash flow over its life:
End of Year 1: $100.00
End of Year 2: $120.00
End of Year 3: $125.00
End of Year 4: $150.00
End of Year 5: $180.00
End of Year 6: $200.00
End of Year 7: $220.00
End of Year 8: $250.00
End of Year 9: $290.00
End of Year 10: $300.00
If we use a discount rate of 10%, what is the present value of the cash flow?
For those unfamiliar with present value, the math behind the calculation is:
PV = -1,000 + 100/(1.1^1) +...

This is July 4th 2009 on the Overlander Trail.
Here's today's route:
Here's T ready to go
The trail is just lined with wood lillies:
And Roses
But, on with the show:
Look the other way:
This is T's river:
And here we are:

In my last post, we wrote a bunch of stock prices for four companies to our database.
We now want to read them back. Also we want to figure out how to use Fluent NHibernate's auto loading of child collections. That is, how to load up a company's prices and the same time we load a company.
Here is the base program I have come up with:
----------------------------------------------
#light
open System
open System.Collections.Generic
open System.IO
open StockPrices
open FluentNHibernate.AutoMap
open FluentNHibernate
let properties = new Dictionary<string, string>()
let connString = "server='BIG_ROCK\LOGGERSEDGE';Initial Catalog=SMDATA;User ID=sa;Password=XXXXX"
properties.Add("connection.provider", "NHibernate.Connection.DriverConnectionProvider")
properties.Add("dialect", "NHibernate.Dialect.MsSql2000Dialect")
properties.Add("connection.driver_class", "NHibernate.Driver.SqlClientDriver")
properties.Add("show_sql", "true")
properties.Add("connection.connection_string", connString)
let autoMappings = (AutoPersistenceModel.MapEntitiesFromAssemblyOf<StockPrices.COMPANY>()).Where(fun t -> if t.Namespace = "StockPrices.StockPrices" then true else false)
let aConfig...

Now that we can read from our database, the next step is to save our price data using FluentNHibernate. The first thing, we need to do is to add a simple table to store our stock price data:
CREATE TABLE [STOCK_PRICE] (
[Id] [int] IDENTITY (1, 1) NOT NULL ,
[COMPANY_ID] [int] NOT NULL ,
[CLOSING_PRICE] [float] NOT NULL ,
[VOLUME] [float] NOT NULL ,
[THE_DATE] [datetime] NOT NULL ,
CONSTRAINT [PK_STOCK_PRICES] PRIMARY KEY CLUSTERED
(
[Id]
) ON [PRIMARY] ,
CONSTRAINT [FK_STOCK_PRICES_COMPANY] FOREIGN KEY
(
[COMPANY_ID]
) REFERENCES [COMPANY] (
[Id]
)
) ON [PRIMARY]
The only interesting thing of note here is the foreign key to the company table. Next we need to add...