Say No to Loops!

One of the reasons people who are used to C, GAUSS, FORTRAN, etc, find Mathematica a bit strange to begin with is that they are used to "procedural" programming styles instead of "functional" programming. So when they start to program in Mathematica and use the loops and nested If[] statements and what-not, they think Mathematica is slow.

This notebook shows just how much you can speed up your computations in Mathematica by using functional programming instead of a loop.  The only problem is you have to be a bit more lateral in your approach to the problem.

The Problem: generating artificial ARMA time series.

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").

Suppose you want to create a function that generates artificial ARMA processes.  You'd probably do something like this:

Version 1: Use a loop.

ARMAlist[length_Integer, sv_, sd_Real, pcoefs_List, qcoefs_List]:=If[Length[Dimensions[pcoefs]] == 1 && 
            Length[Dimensions[qcoefs]] == 1,Module[\({shocklist, maparts, q, 
				p, qplus, thelist}\),
			p=Length[pcoefs];
			q=Length[qcoefs];
			qplus= Reverse[Join[{1}, qcoefs]];
			shocklist=RandomArray[NormalDistribution[0, sd], 
                    length + Max[p, q] ];
			maparts=\(Expand/@\ 
					
                    Table[qplus.shocklist
                          [[Range[t - q, t]
                            ]], {t, q + 1, length + q}]
                      \);
			thelist=
				Join[Table[sv, {p}], Table[0, {length}]];\(Do[thelist = 
                    ReplacePart[thelist, 
			
                      pcoefs.Reverse[
                            thelist
                              [[Range[i - p, i - 1]
                                ]]] + 
                        maparts
                          [[i - p]], 
                      i], {i, p + 1, length + p}]\);Drop[thelist, p] ],Message[ARMAlist::"badcf"]]

Here are the timings for generating 1000 to 12 000 artificial data points on a Power Mac G3 (233 MHz):

timingtesttable1 = 
      Table[Timing[
          ARMAlist[1000 i, 0, 0.5, {0.9, 0.09}, {0.8, 0.8, 0.5}]; ], {i, 1, 
          8}]

ListPlot[(First/@timingtesttable1)/Second]

Ok, so this first version looks OK, runs reasonably well.  But there is a better way.  Functional programming is quicker than that loop.

Version 2: Functional Programming

\(ARMAList[n_Integer, sv_, sd_Real, pcoefs : {_, __}, qcoefs_?VectorQ] /; 
        \ 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}]&, 
                  Table[sv, {p}], mabit])\)]\ ]]\)

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]

There!  Heaps faster and approximately linear in the number of data points to be generated.  The previous version would have taken AGES to generate 10,000 data points, especially on a slower machine, but with the functional programming style you can do it in less than two seconds.  Consider that the 8000-point case is nearly TWENTY times faster!!  And there are fewer lines of code. All right!

One final thing: you might have noticed that I've defined pcoefs in the new version to match {_,__}, that is, a vector with two or more elements.  That's because the AR(1) case admits a simpler function and I decided to code that as a special case.

Luci Ellis
8 July 1998