TrenchInv.m Package Documentation

Installation

The files trenchM.exe and trenchQFM.exe should be copied to the system Mathematica directory where Mathematica.exe is located.  On my computer this path is:

C:\Program Files\Wolfram Research\Mathematica\5.2

The file TrenchInverse.m should be copied to the Mathematica Applications directory.  Typically the path for this is:

C:\Program Files\Wolfram Research\Mathematica\5.2\AddOns\Applications

Package Overview

({{trenchM, \(C\ interface\ function\  - \ \(\(not\)\(\ \)\(normally\)\(\ \)\(used\)\(\ \)\(directly\)\(\ \)\)\)}, {trenchQFM, C interface function  -  
              not normally used directly}, {TrenchInverse, inverse of a p . d .  Toeplitz matrix}, {TrenchLoglikelihood, loglikelihood using Trench algorithm}, {DLLoglikelihood, direct computation of loglikelihood}, {TrenchMean, MLE for μ}, {TrenchInnovationVariance, MLE for (Subscript[σ, z]^2)}, {toe, form Toeplitz matrix from vector}, {toeQ, test is it is a Toeplitz matrix}})

The package MathLink to run the C code programs trenchM.exe and trenchQFM.exe.

The Mathematica function TrenchInverse is used for computing the matrix inverse of a positive-definite Toeplitz matrix.  The package uses MathLink to run a C program to make the computation as efficient as possible.  TrenchInverse requires O((n^2)) flops for an order n matrix and is much faster than Mathematica's Inverse function especially when n is large.

The Mathematica function TrenchLoglikelihood is used for computing the exact loglikelihood.  This function using the MathLink interface to run trenchQMF.exe which uses the Trench algorithm to evaluate the quadratic form (z^′) (Γ^(-1)) z and determinant |Γ|,  where z is a vector of length n and Γ is a positive-definite Toeplitz matrix.   One application of trenchQMF.exe is to the evaluation of the likelihood function of an arbitrary stationary process.  Denoting the autocovariance of this process by (Subscript[γ, k]),  k = 0, 1,  ...  then the covariance matrix of n successive observations is given by \(TraditionalForm\`Γ = \(\(\((γ\_\(i - j\))\)\_\(n×n\)\)\(.\)\)\)  After maximizing over the innovation variance parameter the concentrated loglikelihood function can be written L = -(1/
                  2) log(|Γ|) - 
          (z^′) (Γ^(-1)) z.  Note that this function is unchanged if we use the autocorrelation function (Subscript[ρ, k]) = (Subscript[γ, k])/(Subscript[γ, 0]) instead.  As noted in our paper, the algorithm trenchQMF.exe assumes the autocorrelation form.

As a check on our program we have also written a separate Mathematica function DLLoglikelihood which also computes the exact concentrated loglikelihood but it uses another approach based on the Durbin-Levinson algorithm to obtain the sequence of one-step prediction errors.  This function is entirely written in Mathematica with no use of C code and is used to validate the output of TrenchLoglikelihood.

TrenchMean computes the exact MLE, \(TraditionalForm\`\(\(μ\&^\)\(,\)\)\) for mean given n successive observations from a stationary time series with autocovariance or autocorrelation function c,

\(TraditionalForm\`μ\&^ = \(\(1\_n\%′\) \(Γ\_n\%\(-1\)\) z\)\/\(\(1\_n\%′\) \(Γ\_n\%\(-1\)\) 1\)\),

where z=((Subscript[z, 1]),  ... , (Subscript[z, n])) isthe data vector, (Subscript[1, n]) is a vector of n 1's and (Subscript[Γ, n]^(-1)) is the inverse of the Toeplitz matrix of the vector ((Subscript[c, 0]),  ... , (Subscript[c, n - 1])).

The Mathematica function TrenchInnovationVariance computes the exact MLE for the innovation variance given n successive observations from a stationary time series with autocovariance function c.  The exact MLE for (Subscript[σ, a]^2) given the MLE for c can be written, \(TraditionalForm\`\(\(\(σ\&^\)\_a\%2\)\(=\)\)\)(z^′) (Γ^(-1)) 
        z/n.

Load package

<< TrenchInverse.m

Matrix Inverse

?TrenchInverse

?toe

n = 100 ;phi1 = 0.8 ;c = 
        phi1^Range[0, n - 1]/(1 - phi1^2) ;  (* 
      autocovariance function of AR (1) with parameter phi1 *)

g  =  toe[c];

gi = TrenchInverse[g];

Max[Abs[Flatten[gi - Inverse[g]]]]

The maximum absolute error is very small so we conclude TrenchInverse is working.

Loglikelihood Computation

SeedRandom[231417]; n = 100; phi1 = 0.6;

z = Table[Random[] - 0.5, {n}];

c = phi1^Range[0, n - 1]/(1 - phi1^2) ;  (* 
      autocovariance function of AR (1) with parameter phi1 *)

?TrenchLoglikelihood

TrenchLoglikelihood[c, z]

?DLLoglikelihood

DLLoglikelihood[c, z]

MLE for the Mean and Innovation Variance

Given z=((Subscript[z, 1]),  ... , (Subscript[z, n])) and the autocorrelation function (Subscript[ρ, k])=(Subscript[γ, i]), k = 0, 1,  ... , n - 1the exact mle for the mean is,

\(TraditionalForm\`μ\&^ = \(\(1\_n\%′\) \(Γ\_n\%\(-1\)\) z\)\/\(\(1\_n\%′\) \(Γ\_n\%\(-1\)\) 1\)\)

Note that rescaling the autocovariance function by dividing by (Subscript[γ, 0]) has no effect on \(TraditionalForm\`\(\(μ\&^\)\(.\)\)\)

Given a mean-corrected series z=((Subscript[z, 1]),  ... , (Subscript[z, n])) and (Subscript[γ, i]), k = 0, 1,  ... , n - 1the exact mle for the innovation variance is,

\(TraditionalForm\`\(σ\&^\)\_a\%2 = \(z\^′\) \(Γ\_n\%\(-1\)\) z/n\)

Note that it is assumed that (Subscript[γ, k]) has been rescaled by dividing by the innovation variance.

The following code simulates an AR(1) time series with parameter ϕ and unit innovation, (Subscript[z, 1])~N(
          0, (Subscript[σ, a]^2)/(1 - (ϕ^2))) and (Subscript[z, t]) = ϕ (Subscript[z, t - 1]) + (Subscript[a, t]) where (Subscript[a, t])~NID(0, (Subscript[σ, a]^2)). In the example below, (Subscript[σ, a]^2) = 9.

<< Statistics`NormalDistribution`

ϕ = 0.5; sigsq = 9.0; g0 = sigsq/(1 - (ϕ^2)); n = 800; 
    c = Table[g0 (ϕ^k), {k, 0, n - 1}];
sd = Sqrt[g0];
sig = Sqrt[sigsq];
SeedRandom[2313317];
z1 = Random[NormalDistribution[0, sd]];
z = 
        FoldList[#1*ϕ + #2 &, z1, 
          RandomArray[NormalDistribution[0, sig], 
            n - 1]
];

?TrenchMean

TrenchMean[c, z]

Mean[z]

?TrenchInnovationVariance

TrenchInnovationVariance[c/sigsq, z]

Toeplitz matrix test

R = toe[Range[3]]

toeQ[R]

TrenchInverse[R]