FGN.m Package Documentation

Installation

The file FGN.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

The package FGN.m.m requires that the TrenchInv.m package be installed.  The principal functions in this package are listed in the table below:

({{  FGNAcf, \(theoretical\ autocorrelation\ function\ of\ FGN\)}, {  FGNLL, loglikelihood of FGN}, {   FitFGN, fits FGN time series}, {   FitRegressionFGN, \(\(fits\)\(\ \)\(multiple\)\(\ \)\(linear\)\(\ \)\(regression\)\(\ \)\(with\)\(\ \)\(FGN\)\(\ \)\(error\)\(\ \ \ \)\)}, {   HurstK, Hurst nonparametric estimate of H}, {   SimulateFGN, \(fast\ simulation\ of\ FGN\ using\ FFT\ method\)}})

Other helper functions of interest include:

({{SimulateGLP, simulates general linear process}, {Acvf, sample autocovariance function}, {Acf, sample autocorrelation function}, {StepIAXMatrix , simple step intervention design matrix}, {ProfileLogLRegressionFGN, profile loglikelihood}})

Fractional Gaussian Noise (FGN) Definition

A discrete time series (Subscript[z, t]),  t = 1,  ... , n may be defined as a Gaussian time series with autocorrelation function \(TraditionalForm\`ρ\_k = \((\(|\)\(k + 
                1\)\( | \^\(2  
                    H\)\)\(\(\(-2\)\ k\^\(2  H\)\)\(+\)\)\(|\)\(k - 
                1\)\( | \^\(2  H\)\))\)/2\).  The FGN model was suggested for use with geophysical time series.  An historical review and the basic properties of this model were discussed by Hipel and McLeod (2005).   The maximum likelihood estimator for H and the innovation variance are easily obtained using the TrenchInv.m  package.

Hurst K Estimator

Let (Subscript[z, 1]),  ... , (Subscript[z, n]) be the observed time series.  Define the cumulative range statistic,

\(S\_k = ∑\+\(t = 1\)\%k\((z\_t - z\&_)\)\),k = 1,  ... , n

Let (Subscript[R, n]) = (max_k) (Subscript[S, k]) - (min_k) (Subscript[S, k]) and

\(TraditionalForm\`c\_0 = ∑\+\(t = 1\)\%n\((z\_t - z\&_)\)\^2\)

Then the rescaled adjusted range is \(TraditionalForm\`\(R\_n\%*\) = R\_n/\@c\_0\) and the Hurst K estimator (Hipel and McLeod, p.332),

\(TraditionalForm\`\(R\_n\%*\) = \((n/2)\)\^K\)

or

\(TraditionalForm\`K = \(log(\(R\_n\%*\))\)/\(log(n/2)\)\)

Simulation

Let (Subscript[X, t]),  t = 0, 1,  ... , n be a stationary time series with autocovariance function (Subscript[γ, k]).  Form the sequence,

c=\((γ\_0,  ... , γ\_n, γ\_\(\(n - 1\)\(,\)\) ... , γ\_1)\).

This sequence has length N = 2  n.  Let (Subscript[g, k]) be the finite Fourier transform of c.  That is,

\(TraditionalForm\`g\_k = ∑\+\(j = 0\)\%\(N - 1\)c\_j\ exp {\(-2\) π\ i\ k\ j/N}\)

Let (Subscript[Z, k]), k = 0,  ... , N - 1 be iid complex normal random variables.  Set (Subscript[Z, 0]) and (Subscript[Z, n]) as mean zero and variance 2.  For (Subscript[Z, k]),  k = 1,  ... , n - 1have real and imaginary parts normal mean zero and variance 1.  Finally set \(TraditionalForm\`Z\_\(2  n - k\) = Z\&_\_k\).  Then set for j = 0,  ... , n

\(TraditionalForm\`X\_j = 
      N\^\(\(-1\)/2\)\ \(∑\+\(k = 0\)\%\(N - 1\)\(Z\_k\) \@g\_k\ \ exp {2  π\ i\ k\ j/N}/\@2\)\)

Then (Subscript[X, 0]), (Subscript[X, 1]),  ... ,  (Subscript[X, n]) has mean zero and theoretical autocovariance function, (Subscript[γ, k]).

SimulateGLP implements this method for a general linear process (GLP).  The functions Acf and Acvf for computing the sample autocorrelation and autocovariances functions are included to provide a check on the output from SimulateGLP.  For large n, the input theoretical autocovariance function should agree with the sample autocovariance of the output.

The function  SimulateFGN uses FGNAcf and the function SimulateGLP to simulate FGN.

Reference

Davies,  R. B. and Harte, D. S. (1987).  Tests for Hurst Effect Biometrika  74  95--101.

Illustrative Usage

Load package

<< FGN.m

Average annual riverflow, Nile River at Aswan, 1870-1945.

Series length is n = 75 and we assume the intervention occured at T = 33.   Hipel and McLeod (2005, §19.2.4) fit a step intervention model with T = 33 and an AR(1) error.

nile = {3958.04, 3369.69, 3485.24, 3437.69, 3702.35, 3817.61, 
          2875.58, 3054.69, 4724.15, 3834.01, 3076.77, 2965.76, 3461.71, 
          3141.01, 3371.24, 2988.43, 3607.54, 2946.08, 2709.2, 3294.85, 
          3556.62, 3653.93, 3846.06, 3713.64, 4252.31, 3657.5, 3639.37, 
          3197.72, 3112.75, 2353.68, 2843.65, 2194.93, 2689.43, 2950.91, 
          2247.88, 2628.28, 2491.13, 2792.63, 3321.47, 3058.06, 2889.85, 
          2495.27, 1648.82, 1981.96, 2411.07, 3035.2, 3556.13, 3261.96, 
          2377.89, 2394.96, 2500. , 2610.24, 2743.63, 2744.12, 2338.64, 
          2494.98, 2474.44, 2446.37, 2963.06, 2732.25, 2205.15, 2681.81, 
          2580.54, 2954.38, 3025.94, 2902.78, 2642.46, 2860.24, 2665.41, 
          2306.91, 1848.09, 2569.54, 2503.95, 2438.75, 2211.13};

Reference

Hipel, K.W. and McLeod, A.I., (2005). Time Series Modelling of Water Resources and Environmental Systems  Electronic reprint of our book orginally published in 1994. http://www.stats.uwo.ca/faculty/aim/1994Book.

Nile Minima, 622 AD to 1284 AD.

The original source for this data is Toussoun (1925).  It was reprinted in Beran (1994) and also used as an example by Percival and Walden (2000, Figure 192).  There are n = 663 observations.  Units are in meters.  Sampling interval is 1 year.

NileMin = {11.57, 10.88, 11.69, 11.69, 9.84, 13.22, 11.78, 11.03, 
          12.11, 12.92, 11.24, 11.71, 11.33, 12.27, 11.42, 12.16, 12.59, 
          12.99, 12.32, 11.17, 11.55, 12.32, 10.83, 10.2, 13.94, 11.96, 
          11.48, 10.83, 11.89, 11.33, 10.34, 11.57, 10.34, 10.97, 12.99, 
          11.57, 11.3, 11.55, 13.49, 12.32, 11.03, 11.03, 10.83, 10.27, 
          11.66, 11.48, 12.5, 11.55, 10.47, 10.54, 10.18, 11.89, 11.26, 12.5, 
          12.97, 11.78, 10.43, 11.03, 12.5, 12.72, 11.69, 10.04, 10.83, 
          11.64, 11.24, 10.27, 9.95, 11.69, 12.7, 10.11, 12.47, 11.01, 10.04, 
          10.04, 10.65, 12.23, 11.84, 12.16, 11.8, 11.42, 12.77, 12.06, 
          10.76, 10.76, 11.89, 11.21, 11.78, 10.31, 10.76, 11.78, 12.09, 
          10.22, 12.2, 10.7, 11.26, 10.58, 12.16, 13.58, 11.84, 10.83, 10.97, 
          11.19, 10.97, 10.97, 11.53, 11.53, 11.51, 11.51, 11.51, 11.84, 
          10.97, 10.43, 10.43, 10.02, 11.52, 10.97, 10.34, 10.02, 9.89, 
          10.92, 11.15, 11.15, 10.47, 10.4, 10.38, 10.85, 11.26, 10.58, 
          10.67, 11.15, 12.63, 11.24, 11.1, 10.97, 10.97, 11.57, 10. , 9.91, 
          9.95, 10.13, 10.07, 9.71, 9.71, 9.8, 9.93, 10.43, 10.97, 9.82, 
          9.71, 9.71, 10.65, 10.22, 10.29, 9.89, 10.29, 9.95, 9.82, 10.9, 
          9.8, 9.71, 9.57, 9.89, 9.66, 9.89, 10.22, 10.74, 11.1, 11.1, 10.61, 
          11.51, 11.28, 10.74, 10.43, 10.34, 10.74, 9.66, 10.27, 10.29, 
          10.34, 10.65, 9.89, 10.34, 10.02, 11.28, 11.78, 10.97, 11.42, 
          14.66, 10.97, 11.37, 10.97, 12.59, 13.13, 11.73, 11.69, 11.73, 
          10.88, 11.91, 11.46, 11.6, 11.42, 11.28, 11.69, 11.62, 11.15, 
          11.64, 10.88, 10.79, 10.83, 10.43, 11.1, 10.92, 11.1, 10.47, 10.76, 
          11.1, 10.43, 11.03, 10.34, 10.74, 10.52, 10.11, 10.97, 10.92, 11.1, 
          11.15, 10.97, 11.96, 11.15, 11.62, 11.51, 11.42, 11.26, 11.08, 
          11.87, 11.91, 11.53, 12.54, 11.87, 11.96, 13.31, 14.12, 13.49, 
          12.9, 12.11, 12.32, 11.66, 11.24, 11.46, 10.79, 11.08, 10.97, 
          11.06, 10.72, 10.65, 11.28, 13.4, 9.59, 9.59, 11.37, 11.33, 11.37, 
          11.51, 11.17, 11.57, 11.57, 11.33, 11.1, 11.55, 11.89, 12.6, 11.89, 
          11.51, 10.97, 12.09, 11.3, 12.95, 13.08, 12.5, 12.05, 13.1, 12.5, 
          11.55, 11.01, 11. , 11.03, 11.21, 11.21, 10.97, 11.06, 12.59, 
          12.61, 11.24, 11.96, 12.05, 12.05, 11.19, 10.88, 12.5, 10.94, 
          11.98, 11.21, 11.64, 12.11, 11.53, 11.46, 11.26, 12.88, 11.75, 
          11.71, 10.81, 11.33, 11.64, 11.55, 11.55, 11.55, 11.6, 10.94, 
          10.54, 10.67, 10.44, 9.48, 10.99, 10.16, 10.65, 10.67, 10.72, 
          10.76, 10.81, 11.96, 11.96, 11.51, 10.88, 11.28, 11.51, 12.36, 
          12.16, 12.88, 12.97, 11.82, 13.06, 10.43, 11.84, 10.54, 11.69, 
          10.43, 9.8, 10.72, 11.89, 11.51, 11.42, 11.93, 11.51, 10.97, 11.44, 
          10.97, 10.94, 11.53, 11.08, 9.35, 10.81, 10.81, 10.97, 11.46, 12.5, 
          11.51, 10.43, 10.43, 10.43, 10.7, 11.24, 11.37, 11.46, 10.99, 
          10.54, 10.45, 10.7, 11.42, 10.74, 11.01, 12.2, 11.96, 10.97, 12.07, 
          11.19, 11.6, 11.51, 10.25, 10.97, 11.37, 10.07, 10.34, 10.43, 
          10.43, 9.8, 10.79, 11.69, 12.5, 13.24, 12.09, 11.42, 10.61, 10. , 
          10.88, 11.28, 11.42, 12.59, 11.42, 11.48, 10.88, 11.42, 11.19, 
          11.3, 10.88, 12.5, 11.37, 11.08, 11.1, 11.73, 11.73, 11.96, 11.89, 
          12. , 13.51, 12.74, 12.27, 13.1, 11.48, 11.51, 11.51, 11.82, 11.82, 
          11.51, 11.33, 11.3, 11.51, 11.66, 10.7, 12. , 10.74, 11.1, 12.92, 
          11.78, 11.28, 10.97, 13.04, 11.03, 12.59, 11.19, 11.19, 11.19, 
          10.81, 11.96, 10.85, 11.01, 11.03, 11.46, 12.11, 11.69, 11.44, 
          11.91, 11.89, 11.82, 12.43, 12.43, 12.27, 11.89, 11.91, 11.55, 
          12.09, 12.18, 12.11, 12.09, 11.64, 11.35, 11.21, 11.37, 12.54, 
          14.57, 12.99, 12.77, 12.77, 11.78, 12.7, 13.13, 13.33, 12.7, 12.45, 
          12.45, 12.11, 12.65, 13.46, 13.46, 12.9, 12.95, 12.86, 12.59, 
          12.54, 14.21, 12.68, 12.63, 13.35, 13.13, 12.65, 13.19, 13.51, 
          12.77, 13.17, 12.68, 12.63, 11.12, 12.07, 12.92, 12.05, 12.23, 
          12.05, 11.53, 11.82, 12.45, 12.05, 11.51, 10.79, 11.51, 10.81, 
          11.28, 12.09, 11.57, 12.77, 12.59, 12.09, 12.2, 11.84, 12.2, 11.93, 
          12.47, 12.52, 12.59, 12.99, 11.73, 11.82, 11.8, 11.8, 13.31, 12.07, 
          12.36, 11.51, 11.82, 11.42, 11.91, 12.59, 11.66, 11.96, 12.41, 
          12.52, 12.41, 12.52, 11.57, 11.26, 11.64, 10.88, 11.73, 12.52, 
          12.88, 13.01, 12.86, 12.23, 12.32, 11.84, 12.07, 12.5, 12.56, 
          12.11, 12.16, 12.09, 12.09, 12.07, 11.51, 10.97, 10.97, 9.89, 9.66, 
          10.47, 10.56, 11.1, 12.9, 11.51, 11.66, 11.96, 11.96, 11.1, 11.1, 
          11.19, 11.19, 10.74, 11.06, 11.28, 12.18, 10.98, 10.44, 10.56, 
          10.58, 10.98, 10.43, 10.38, 11.42, 11.42, 11.93, 11.03, 9.89, 9.36, 
          11.42, 11.42, 11.51, 11.51, 11.8, 12.59, 11.96, 11.42, 11.69, 
          11.96, 11.42, 11.28, 10.43, 10.97, 11.42, 12.05, 12.05, 11.64, 
          11.6, 11.96, 11.12, 11.69, 11.1, 11.78, 11.33, 11.53, 11.39, 11.55, 
          11.87, 11.96, 12.2, 11.66, 11.28, 11.01, 11.57, 11.75, 11.42, 
          11.87, 12.54, 11.98, 12.63, 12.83, 12.52, 11.6, 12.34, 12.34, 
          12.32, 13.06, 12.05, 10.54, 11.51, 11.08, 10.97};

Reference

Beran, J. (1994).  Statistics for long-memory processes.  Chapman & Hall, New York.

Percival and Walden (2000).  Wavelet methods for time series analysis.  Cambridge University Press.

Autocorrelation Function

?FGNAcf

Table[FGNAcf[k, 0.7], {k, 1, 20}]

Loglikelihood

?FGNLL

y = Table[Random[] - 0.5, {100}];

FGNLL[0.3, y]

Fitting FGN Model

?FitFGN

FitFGN[y]

FitFGN[y, MeanZeroQ -> False]

FitFGN[nile, MeanZeroQ -> False]

Regression With FGN Error

Here we fit an intervention analysis model, (Subscript[z, t]) = μ + ω (Subscript[S, t]^((T))) + (Subscript[e, t]),  
      where μ and ω are the regression parameters, (Subscript[e, t]) is assumed to follow be a FGN process with parameter H and (Subscript[S, t]^((T))) is the unit step function,

(Subscript[S, t]^((T)))= {({{0, t < T}, {1, t ≥ T}})

We assume here T = 33.

?StepIAXMatrix

X = StepIAXMatrix[Length[nile], 33];

X

FitRegressionFGN[X, nile]

Contour Plot of 95% C.I. Using the Likelihood-ratio Method

ansFGN = FitRegressionFGN[ X,  nile];

alphaHat = RegressionCoefficients /. ansFGN

𝔏max = 
      ProfileLogLRegressionFGN[alphaHat, X, nile]

𝔏[μ_, ω_] := 
        Exp[ ProfileLogLRegressionFGN[{μ, ω}, X, 
              nile]  -  𝔏max ];

\(\(gFGN = 
        ContourPlot[𝔏[μ, ω], {μ, 2200, 
            4100}, {ω, 100, \(-1000\)}, Contours -> {0.1465}, 
          ContourShading -> False, 
	PlotPoints -> 50, 
          PlotRange -> All, 
          Epilog -> {AbsolutePointSize[3], RGBColor[0, 0, 0], 
              Point[alphaHat]}, \ 
          FrameLabel -> {"μ", "ω"}, 
          RotateLabel -> False, 
          ContourStyle -> {RGBColor[0, 1, 0]}];\)\)

SimulateGLP

?SimulateGLP

We simulate an AR(1) with autocovariance function (Subscript[γ, k]) = (ϕ^k)/(1 - (ϕ^2)),  
      k ≥  0 and compare (Subscript[γ, k]) with the sample autocovariance function of the output,

n = (2^14); phi = -0.8; g0 = 1/(1 - (phi^2)); 
    ac = Table[g0* (phi^k), {k, 0, n}];

z = SimulateGLP[ac];

Short[z, 5]

Autocovariance Functions

?Acvf

ta = Take[ac, 11];

sa = Acvf[z, MaxLag -> 10]

\(TableForm[{Range[0, 10], ta, sa} // Transpose, 
      TableHeadings -> {None, {"Lag", "Theoretical", "Sample"}}]\)

SimulateFGN

We simulate a FGN with autocovariance function equal to the FGN autocorrelation function and compare (Subscript[γ, k]) with the sample autocovariance function of the output,

?SimulateFGN

z = SimulateFGN[0.75, (2^14)];

Short[z, 5]

Autocovariance Functions

ta = Table[FGNAcf[k, 0.75], {k, 0, 10}]

sa = Acvf[z, MaxLag -> 10];

\(TableForm[{Range[0, 10], ta, sa} // Transpose, 
      TableHeadings -> {None, {"Lag", "Theoretical", "Sample"}}]\)

Further Simulation Examples

Simulate and estimate H using the Hurst K method. Calculate mse

\(H0 = 0.9; 
    H = Table[z = SimulateFGN[H0, 1000]; 
HurstK[z], {10}]\)

Mean[(((H - H0))^2)]

Simulate and estimate H exact MLE

\(H0 = 0.5; 
    H = Table[z = SimulateFGN[H0, 300]; 

        HurstCoefficient /. Rest[FitFGN[z]], {10}]\)

Mean[(((H - H0))^2)]

\(H0 = 0.7; 
    H = Table[z = SimulateFGN[H0, 300]; 

        HurstCoefficient /. Rest[FitFGN[z]], {10}]\)

Mean[(((H - H0))^2)]