Generating time series in Mathematica

First, load in the required packages

It's important because otherwise you can't generate normally distributed shocks

Needs["Statistics`ContinuousDistributions`"]

This one is also useful in general, and the function EGARCHList below depends on some of its contents.

Needs["LinearAlgebra`MatrixManipulation`"]

Note that this package also loads in some others, as shown below.  This is important if you want to convert this notebook into package form.

$Packages

Random Walk —with or without drift

This is a very simple function, using the functional programming construct NestList[] with a simple pure function. Note that the length of the series generated is constrained to be an integer, while the variables drift and sd -- the drift and the standard deviation of the shocks -- have default values.  I haven't constrained the starting value (sv) to be numerical -- sometimes you might want to explore the properties of something in algebraic terms.

randomwalk[n_Integer, sv_, sd_:  1, drift_:  0] := 
      Rest[NestList[#1 + Random[NormalDistribution[0, sd]] + drift &, sv, n]]

ARMA: generate artificial time series of any order

Anyone who's done econometrics or time series work would be familiar with this sort of process:

	\(TraditionalForm
      \`y\_t = ∑\+\(i = 1\)\%p\( ϕ\_i\) y\_\(t - 1\)\  + 
          ∑\+\(j = 0\)\%q\( θ\_j\) ϵ\_\(t - j\)\)

The ϵ's are usually normally distributed, or at least iid. (should that be "normally normally").

\(ARMAList[n_Integer, sv_, sd : \((_Real | _Integer)\), pcoefs : {_, __}, 
          qcoefs_?VectorQ, const_:  0] /; \ VectorQ[pcoefs] := \ 
      With[{p = Length[pcoefs], q = Length[qcoefs], 
          macoefs = Flatten[Reverse[Join[{1}, qcoefs]]], 
          rpcoefs = Reverse[pcoefs], 
          rawshocks = 
            RandomArray[NormalDistribution[0, sd], n + Length[qcoefs]]}, 
        With[{mabit = 
              Expand\ /@
                \((Take[
                    \((macoefs.
                        \((\(RotateLeft[rawshocks, #]&\)\ /@\ Range[0, q])\))
                      \), n])\)}, 
			
          Rest[Last\ /@\ 
              \((FoldList[Join[Rest[#1], {rpcoefs.#1\  + #2 + const}]&, 
                  Table[sv, {p}], mabit])\)]\ ]]\)

Special Cases for AR(1) processes.  The pattern-matching engine will pick this up automatically, so you can put the single AR coefficient in list brackets or not as you see fit.

\(ARMAList[n_Integer, sv_, sd : \((_Real | _Integer)\), pcoefs : {_}, 
          qcoefs_?VectorQ, const_:  0] /; \ VectorQ[pcoefs] := \ 

      ARMAList[n, sv, sd, pcoefs[[1]], 
        qcoefs, const]\)

\(ARMAList[n_Integer, sv_, sd : \((_Real | _Integer)\), 
        pcoef : \((_Real | _Integer | \ _Symbol)\), qcoefs_?VectorQ, 
        const_:  0] := \ 
      With[{q = Length[qcoefs], 
          macoefs = Flatten[Reverse[Join[{1}, qcoefs]]], 
          rawshocks = 
            RandomArray[NormalDistribution[0, sd], n + Length[qcoefs]]}, 
        With[{mabit = 
              Expand\ /@
                \((Take[
                    \((macoefs.
                        \((\(RotateLeft[rawshocks, #]&\)\ /@\ Range[0, q])\))
                      \), n])\)}, 
			
          FoldList[pcoef\ *\ #1\  + #2 + const\ &, sv, mabit]\ ]]\)

Special Cases for shock vectors that are already supplied.

\(ARMAList[ϵ_?VectorQ, sv_, pcoefs : {_, __}, qcoefs_?VectorQ, 
          const_:  0] /; \ VectorQ[pcoefs] := \ 
      With[{p = Length[pcoefs], q = Length[qcoefs], 
          macoefs = Flatten[Reverse[Join[{1}, qcoefs]]], 
          rpcoefs = Reverse[pcoefs], 
          T = Length[ϵ] - Length[qcoefs]}, 
        With[{mabit = 
              Expand\ /@
                \((Take[
                    \((macoefs.
                        \((\(RotateLeft[ϵ, #]&\)\ /@\ Range[0, q])\))
                      \), T])\)}, 
			
          Rest[Last\ /@\ 
              \((FoldList[Join[Rest[#1], {rpcoefs.#1\  + #2 + const}]&, 
                  Table[sv, {p}], mabit])\)]\ ]]\)

\(ARMAList[ϵ_?VectorQ, sv_, pcoef : {_}, qcoefs_?VectorQ, 
          const_:  0] /; \ VectorQ[pcoefs] := \ 

      ARMAList[ϵ, sv, 
        pcoef[[1]], qcoefs, const]\)

\(ARMAList[ϵ_?VectorQ, sv_, 
        pcoef : \((_Real | _Integer | \ _Symbol)\), qcoefs_?VectorQ, 
        const_:  0] := \ 
      With[{q = Length[qcoefs], 
          macoefs = Flatten[Reverse[Join[{1}, qcoefs]]], 
          T = Length[ϵ] - Length[qcoefs]}, 
        With[{mabit = 
              Expand\ /@
                \((Take[
                    \((macoefs.
                        \((\(RotateLeft[ϵ, #]&\)\ /@\ Range[0, q])\))
                      \), T])\)}, 
				
          FoldList[pcoef\ *\ #1\  + #2\  + const\ &, sv, mabit]\ ]]\)

Here is an example.

ListPlot[ARMAList[1000, 0, 0.5, {0.9, 0.09}, {0.8, 0.8, 0.5}, 0.2], 
        PlotJoined -> True, Frame  ->  True, 
        PlotStyle -> {AbsoluteThickness[2], Hue[0.57]} ] // Timing

Amazingly, this function is linear in the length of the series to be created.  Not many other functions manage that!

timingtesttable = 
      Table[Timing[
          ARMAList[1000 i, 0, 0.5, {0.9, 0.09}, {0.8, 0.8, 0.5}]; ], {i, 1, 
          30}]

ListPlot[(First/@timingtesttable)/Second]

ARCH: generate artificial time series with autoregressive conditional heteroskedasticity

These functions are based on the time series generation procedures written by Pedro J F de Lima of the Johns Hopkins University in the GAUSS programming language, available at ftp://speedy.econ.jhu.edu/public/gauss/st_proc.prc.  Note that the function tests to ensure that the coefficients are standardised, ie that ω =1 – ∑α – ∑β.  If this is not true the function definition will not match and nothing will happen.

With both AR and MA components, we have GARCH:

\(ARCHList[n_Integer, ω_Real, α : {__Real}, 
          β : {__Real}] /; \ 
        \((ω == 1 - \((Plus@@α)\) - \((Plus@@β)\))\) := 

      	With[{m = Max[Length[α], Length[β]]}, 
	
        With[{q = Length[α], p = Length[β], 
            arev = Reverse[α], brev = Reverse[β], 
            z = RandomArray[NormalDistribution[0, 1], m + n]}, 
			
          Module[{h = Table[1, {n + m}], 
              e = Join[z[[Range[m]]], 
                  Table[0, {n}]]}, 
				
            Do[h[[i]] = 
                ω\  + 
                  \((\((e[[Range[i - q, i - 1]
                              ]]\ )\)^2)\).arev + 
                  h[[Range[i - p, i - 1]
                        ]]\ .brev; 
				
              \(e[[i]] = 
                z[[i]] + 
                  Sqrt[h[[i]]]; 
              \), {i, m + 1, n + m}]; 
				{Drop[h, m], Drop[e, m]\ }]]]
        \)

Without either AR or MA components we have white noise.

ARCHList[n_Integer, ω_Real, {}, {}] := ARCHList[n, ω]

ARCHList[n_Integer, ω_Real] := {
        RandomArray[NormalDistribution[0, 1], n], Table[1, {n}]}

With only AR components, it's ARCH but not GARCH.

ARCHList[n_Integer, ω_Real, α : {__Real}, {}] := 
      ARCHList[n, ω, α]

\(ARCHList[n_Integer, ω_Real, α : {__Real}] /; 
        \((ω == 1 - \((Plus@@α)\))\) := 
		
      With[{q = Length[α], arev = Reverse[α], 
          z = RandomArray[NormalDistribution[0, 1], q + n]}, 
			
        Module[{h = Table[1, {n + q}], e = z}, 
				
          Do[h[[i]] = 
              ω\  + 
                \((\((e[[Range[i - q, i - 1]
                            ]]\ )\)^2)\).arev; 
				
            \(e[[i]] = 
              z[[i]] + 
                Sqrt[h[[i]]]; 
            \), {i, q + 1, n + q}]; 
				{Drop[h, q], Drop[e, q]\ }]]\)

The following cells generates plot of some GARCH data.

ListPlot[
      ARCHList[1000, 0.05, {0.3, 0.2, 0.15}, {0.2, 0.1}]
          [[1]], PlotJoined -> True, 
      Frame -> True, PlotRange -> All, 
      PlotStyle -> {AbsoluteThickness[2], GrayLevel[0.5]}]

ListPlot[
      ARCHList[1000, 0.05, {0.9, 0., 0.04}, {0.0, 0.01}]
          [[1]], PlotJoined -> True, 
      Frame -> True, PlotRange -> All, 
      PlotStyle -> {AbsoluteThickness[2], GrayLevel[0.5]}]

EGARCH: artificial time series

These functions are based on the time series generation procedures written by Pedro J F de Lima of the Johns Hopkins University in the GAUSS programming language, available at ftp://speedy.econ.jhu.edu/public/gauss/st_proc.prc.  Note that  I have used some of the functions in the standard Mathematica package "LinearAlgebra`MatrixManipulation`".

Needs["LinearAlgebra`MatrixManipulation`"]

$Packages

\(EGARCHList[z : {__Real}, uncmean_, arpar_?VectorQ, 
        mapar_?VectorQ, {θ1_, θ2_}] := 
	
      With[{p = Length[arpar], q = Length[mapar] + 1, 
          const = uncmean + \((1 - Plus@@arpar)\), n = Length[z]}, \ 
		
        Module[{zlag, gz, σ2}, 
			
          zlag = \(Take[RotateLeft[z, #], n]&\)\ /@\ Range[0, q - 1]; 
			
          gz = AppendRows[ZeroMatrix[q], 
              θ1*zlag + θ2\ *\ \((Abs[zlag] - Sqrt[2./Pi])\)]; 
          
			σ2 = 
            ARMAList[const + Reverse[Join[{1}, mapar]].gz, uncmean, 
              arpar, {0}]; 

          			{z*Sqrt[Exp[Drop[σ2, q - 1]\ ]], 
            Drop[σ2, q - 1]}\ ]]\)

Fractional Noise: three methods.

These functions are based on the time series generation procedures written by Pedro J F de Lima of the Johns Hopkins University in the GAUSS programming language, available at ftp://speedy.econ.jhu.edu/public/gauss/gen_fnoi.prc.  There are three methods, finite Fourier transforms of Davies and Harte (1987), Hosking's recursive method (1984) and a method based on the Cholesky Decomposition. ( Note that this requires the standard Cholesky package to be read in first.) The choice between the three methods is handled as an option to the main function, then three subsidiary functions are used.  These implementations are preliminary, as I cannot give any guarantees that I have understood the GAUSS code properly.  In particular, the Hoskings method gives results that look rather different to the other two.

This reads in the package.

Needs["LinearAlgebra`Cholesky`"]

To do the Cholesky decomposition method, you also need this subsidiary function.   
Toeplitz matrices are n × n symmetric matrices with elements (Subscript[a, i, j]) = (Subscript[a, i - 1, j - 1]) for i,j = 2…n.

createToeplitzMatrix[vec_?VectorQ] :=  
      With[{n = Length[vec]}, 
        Array[vec[[Abs[#1 - #2] + 1]]&, {
            n, n}]]

createToeplitzMatrix[{α, β, γ, δ, 
          ϵ}] // MatrixForm

It isn't the fastest, but it works.

Table[
      Timing[createToeplitzMatrix[Table[Random[], {i}]]; ], {i, 50, 300, 50}]

This sets up the FFT method as the default.

Options[generateFractionalNoise] = {Method -> FiniteFourier}

\(generateFractionalNoise[n_Integer, d_Real, opts___?OptionQ] /; 
        \((d > \(-0.5\)\  && \ d < 0.5)\) := 
	
      Module[{way}, 
        way = \(Method /. {opts}\) /. Options[generateFractionalNoise]; 
		
        Switch[way, FiniteFourier, fourierFNoise[n, d], \ Hoskings, 
          hoskingsFNoise[n, d], CholeskiDecomp, cholFNoise[n, d]]\ ]\)

\(fourierFNoise[n_Integer, d_Real] /; 
        \((d > \(-0.5\)\  && \ d < 0.5)\) := 
	
      Module[{gk, u, v, z, x1, auto}, 
        With[{𝓃 = n + 1}, 
          With[{rhopart = 
                Thread[\((Range[𝓃] - 1 + d)\)/
                    \((Range[𝓃] - d)\)], 
              auto0 = Gamma[1 - 2\ d]/\((Gamma[1 - d]^2)\)}, 
            With[{ρ = FoldList[Plus, First[rhopart], Rest[rhopart]]}, 

              				auto = auto0*ρ; 
				
              gk = Re[Fourier[
                    Join[{auto0}, Drop[auto, \(-1\)], 
                      Reverse[Drop[auto, \(-2\)]]]]\ ]; 
				
              u = Thread[
                  Join[{\@2.}, Table[1, {𝓃 - 2}], {\@2.}]*
                    RandomArray[NormalDistribution[0, 1], 𝓃]]; 

              				v = 
                Join[{0.}, 
                  RandomArray[NormalDistribution[0, 1], 𝓃 - 2], {
                    0.}]; 
				
              z = Join[Thread[u + \((I*v)\)], 
                  Reverse[
                    Thread[Drop[Rest[u], \(-1\)] + 
                        \((I*Drop[Rest[v], \(-1\)])\)]]]; 
				
              x1 = 0.5*Re[Fourier[Thread[\@gk*z]]]/\@\(𝓃 - 1\); 

              				Take[x1, 𝓃 - 1] - 
                \((Plus@@Take[x1, 𝓃 - 1])\)/\((𝓃 - 1)\)]]]]
        \)

\(hoskingsFNoise[n_Integer, d_Real] /; 
        \((d > \(-0.5\)\  && \ d < 0.5)\) := 
	
      Module[{vlist, ϕtt, ϕlists}, 
	
        With[{k = Range[n], v0 = Gamma[1 - 2\ d]/\((Gamma[1 - d]^2)\)}, 

          			ϕtt = d/\((k - d)\); 
			
          vlist = Drop[
              FoldList[\((1 - \((#2)\)^2\ )\)*\ #1\ &\ , v0, ϕtt], 
              \(-1\)]; 
			
          ϕlists = 
            FoldList[Flatten[Join[{#1} - #2*\ Reverse[{#1}], {#2}]]&\ , \ 
              First[ϕtt], Rest[ϕtt]]; 
		\ \ 
          Rest[Flatten[
              Last[
						
                FoldList[
                  Join[{#1}, 
                      Flatten[{
                          \@\((#2[[1]])
                                  \)\ *Random[NormalDistribution[0, 1]] + 
                            \((Flatten[{#2[[2
                                      ]]}].Reverse[
                                  Flatten[{#1}]])\)}]]\ &, \ 
                  \@v0*Random[NormalDistribution[0, 1]], 
                  Transpose[{vlist, ϕlists}]]\ ]]]\ ]]\)

\(\(cholFNoise[n_Integer, d_Real] /; \((d > 0.\  && \ d < 0.5)\) := 
	
      With[{τ = Range[n], σe = 1}, 
	
        With[{	auto1 = 
              Exp[LogGamma[1] + LogGamma[1 - 2\ d] - LogGamma[d] - 
                  LogGamma[1 - d] + LogGamma[d + τ] - 
                  LogGamma[1 - d + τ]]}, 
		
          Transpose[
              CholeskyDecomposition[createToeplitzMatrix[auto1]]].RandomArray[
              NormalDistribution[0, 1], n]\ ]]
		\)\)

The Fourier method seems to be fastest.  The Timings below are from a Power Macintosh G3 233 MHz with 64 Mb of RAM and virtual memory on.

ListPlot[generateFractionalNoise[500, 0.49999], PlotJoined -> True]
         // Timing

{0.333333 Second}

ListPlot[generateFractionalNoise[100, 0.49, Method -> Hoskings], 
        PlotJoined -> True] // Timing

{0.233333 Second}

ListPlot[
        generateFractionalNoise[100, 0.49, Method -> CholeskiDecomp], 
        PlotJoined -> True] // Timing

{2.13333 Second}