Cholesky Decomposition

Bayesian Methods for Empirical Macroeconomics
A Course Prepared the Global Initiative Program, Bank of Korea Academy
http://personal.strath.ac.uk/gary.koop/BoK_course.html

#r @"../packages/MathNet.Numerics.3.2.3/lib/Net40/MathNet.Numerics.dll"
#r @"../packages/FSPowerPack.Core.Community.3.0.0.0/Lib/Net40/FSharp.PowerPack.dll"

open System
open System.Collections.Generic
open System.Data
open MathNet.Numerics.Distributions

// http://www.quantstart.com/articles/Cholesky-Decomposition-in-Python-and-NumPy
// http://rosettacode.org/wiki/Cholesky_decomposition#OCaml
// http://www.math.sjsu.edu/~foster/m143m/cholesky.pdf
// http://stackoverflow.com/questions/8565489/improvement-of-matrix-calculation-in-f

// Algorithm
let cholesky (matA:matrix) =
   let x = matA.NumRows
   let mutable s = 0.0
   let res = Matrix.create x x 0.0
   for i in 0 .. x - 1 do
    for j in 0 .. i do
        s <- 0.0  
        for k in 0 .. j - 1 do
            s <- s + (res.[i,k] * res.[j,k])
        if i = j then res.[i,j] <- sqrt(matA.[i,i] - s)
        else res.[i,j] <- ((matA.[i, j] - s) / res.[j,j])  
   res
   
 
//let a =  matrix [  [25.0; 15.0; -5.0];
//                    [15.0; 18.0;  0.0];
//                    [-5.0;  0.0; 11.0] ];

let a = matrix [    [16.0; 4.0;  4.0;  -4.0];
                    [4.0; 10.0;  4.0;  2.0];
                    [4.0; 4.0; 6.0; -2.0];
                    [-4.0; 2.0; -2.0; 4.0] ];

let abc = cholesky a;;

Numerical Integration (Monte Carlo vs Trapezoidal’s Rule)

This is the exercise 2 of question 1, the details please refer the web site of Professor Gary Koop. I found that I have done something wrong, the answer of Monte Carlo is different with Trapezoidal’s Rule. If anyone else find the solution, please feel free to leave me a message.

Bayesian Methods for Empirical Macroeconomics
A Course Prepared the Global Initiative Program, Bank of Korea Academy
http://personal.strath.ac.uk/gary.koop/BoK_course.html

(* ---------------------- *)
(* Question 1 Ex2  *)
(* -------------------------- *)
#r "System.Windows.Forms.DataVisualization.dll"
#r @"../packages/MathNet.Numerics.3.2.1/lib/Net40/MathNet.Numerics.dll"

open System
open System.Net
open System.Windows.Forms
open System.Windows.Forms.DataVisualization.Charting
open Microsoft.FSharp.Control.WebExtensions
open MathNet.Numerics.Distributions
open MathNet.Numerics.Statistics

let theta_mean = 1.0
let theta_var = 4.0

(* the number of draws we obtain will affect the precision *)
(* of the approximation using Monte Carlo integration *)
let number_of_draws = 1000000

//f : function
//a : lower limit
//b : upper limit
//n : count of intervals  
let Integrate (f : float -> float)  (a:float) (b:float) (n:int) =
  let h = (b - a) / float n //step width
  let TrapezoidalRule (f: float -> float) (tp:float*float) = 
    ((snd tp) - (fst tp)) * ((f (fst tp) + f (snd tp)) / 2.0)
  let startPoints = [a..h..b-h]
  let endPoints = [a+h..h..b]
  let intervals = List.zip startPoints endPoints
  List.map (TrapezoidalRule f) intervals |> List.sum

// (x.^2).*(0.1995*exp((-(x-1).^2)/8))
let Numerical_integration = Integrate (fun x -> (x*x)*(0.1995*exp((-((x-1.0)*(x-1.0)))/8.0))) -10.0 10.0 number_of_draws;;

(* ---------------------------------- *)
(* mean value theorem for integration *)
(* ---------------------------------- *)
// A normally distributed random generator
let normd = new Normal(0.0, 1.0)

let theta_post = seq { for i in 1 .. number_of_draws do yield theta_mean + (sqrt(theta_var)* normd.Sample()) }

// Lets estimate the expected value of the function theta^2 using MC integration
let theta_sq = theta_post |> Seq.map (fun x -> x * x)

// Seq.length theta_sq |> printfn "Length: %d"
let MC_integration = theta_sq.Mean()

printf "Numerical integration (Trapezoidal's Rule)"
Numerical_integration;;
printf "     "
printf "Monte Carlo integration"
MC_integration;;
// printf "Monte Carlo integration (alternative method)"
// theta_sq1;;

Randomize number generator in F# (Bayesian computing in functional language. )

Over the last few days, I study econometric with a book “Gary Koop Bayesian Econometrics”. This book is NOT over theoretical/mathematical orientated. It is important for me, there are many book which are related with “Bayesian Econometrics”, but all of them are over mathematical. This book is an exceptional, Professor Gary provided all source code on his web-site. He have provided some question sheets which take reader through some basic programming skills relating to course materials. Anyone else would know more about Bayesian Econometrics, you should take a look on his web-site. You would find/learn more than you expected.

I wrote a F# versions for his source code. The purpose of this article was to show the algorithms, and to borrow an implementation in F# and the examples may be interesting for those who like bayesian computing in functional language.

Professor Gary Koop conducted a course for various central-bank all over the world.
====================================================================================
Bayesian Methods for Empirical Macroeconomics
A Course Prepared the Global Initiative Program, Bank of Korea Academy
http://personal.strath.ac.uk/gary.koop/BoK_course.html

(* ---------------------- *)
(* Question 1 Ex1  *)
(* -------------------------- *)
#r "System.Windows.Forms.DataVisualization.dll"
#r @"../packages/MathNet.Numerics.3.2.1/lib/Net40/MathNet.Numerics.dll"

open System
open System.Net
open System.Windows.Forms
open System.Windows.Forms.DataVisualization.Charting
open Microsoft.FSharp.Control.WebExtensions
open MathNet.Numerics.Distributions
open MathNet.Numerics.Statistics

let num_draws = 10000

(* ---------------------- *)
(* Do results for uniform random number generator *)
(* -------------------------- *)
let uniformRandom =
    printf "RESULTS - MEAN, STD "
    printf "-----------------------------------------------------------------"
    printf " "
    printf "UNIFORM"
    let normd = new Normal(0.0, 1.0)
    let tempp = normd.Samples() |> Seq.take num_draws |> Seq.toList
    let unifans = [tempp.Mean(), tempp.StandardDeviation()]
    let trueans = [0.5, (1.0/sqrt(12.0))]
    printf "%f %f\n",unifans, trueans
    
(* ---------------------- *)
(* Do results for Normal random number generator  *)
(* ---------------------- *)
let standardNormalRandom =
    printf "STANDARD NORMAL"
    let normd = new Normal(0.0, 1.0)
    let tempp = normd.Samples() |> Seq.take num_draws |> Seq.toList
    let normans = [tempp.Mean(), tempp.StandardDeviation()]
    let trueans = [0.0, 1.0]
    printf "%f %f\n",normans, trueans

(* ---------------------- *)
(* Do results for Student-t random number generator  *)
(* ---------------------- *)
let studentTRandom =
    printf "STUDENT-T(3)"
    let stdT = new StudentT(3.0,(float num_draws),1.0)
    let tempp = stdT.Samples() |> Seq.take num_draws |> Seq.toList
    let tans = [tempp.Mean(), tempp.StandardDeviation()]
    let trueans = [0.0, sqrt(3.0)]
    printf "%f %f\n",tans, trueans

(* ---------------------- *)
(* Do results for Beta random number generator *)
(* ---------------------- *)
let BETARandom =
    printf "BETA(3,2)"
    let beta = new Beta(3.0,2.0)
    let tempp = beta.Samples() |> Seq.take num_draws |> Seq.toList
    let betaans = [tempp.Mean(), tempp.StandardDeviation()]
    let trueans = [(3/(3+2)),( 1 / 5 )]
    printf "%f %f\n",betaans, trueans

(* ---------------------- *)
(* Do results for Exponential random number generator *)
(* ---------------------- *)
let exponentialRandom =
    printf "Exponential (5)"
    let exp = new Exponential(5.0)
    let tempp = exp.Samples() |> Seq.take num_draws |> Seq.toList
    let expans = [tempp.Mean(), tempp.StandardDeviation()]
    let trueans = [5.0, 5.0]
    printf "%f %f\n",expans, trueans

(* ---------------------- *)
(* Do results for Chi-Square random number generator *)
(* ---------------------- *)
let chiSquareRandom =
    printf "ChiSquared(3)"
    let chi2 = new ChiSquared(3.0)
    let tempp = chi2.Samples() |> Seq.take num_draws |> Seq.toList
    let chi2ans = [tempp.Mean(), tempp.StandardDeviation()]
    let trueans = [3.0, sqrt(6.0)]
    printf "%f %f\n",chi2ans, trueans

(* ---------------------- *)
(* Do results for Gamma random number generator *)
(* ---------------------- *)
let gammaRandom =
    printf "Gamma(4,2)"
    let gamma = new Gamma(4.0,2.0)
    let tempp = gamma.Samples() |> Seq.take num_draws |> Seq.toList
    let gamans = [tempp.Mean(), tempp.StandardDeviation()]
    let trueans = [8.0, 4.0]
    printf "%f %f\n",gamans, trueans

uniformRandom;;
standardNormalRandom;;
studentTRandom;;
BETARandom;;
exponentialRandom;;
chiSquareRandom;;
gammaRandom;;