Time-series decomposition in Mathematica

Question

I'm studying time-series in R with this book, and there is a nice command in R that creates decompositions. Inside Mathematica 9 the command can be executed as:

Needs["RLink`"]
InstallR[];
REvaluate["{
    url <- \"http://www.massey.ac.nz/~pscowper/ts/cbe.dat\"
    CBE <- read.table(url, header = T)
    Elec.ts <- ts(CBE[, 3], start = 1958, freq = 12)
    plot(decompose(Elec.ts, type = \"multi\"))
 }"]

and creates this plot:

The first part of my attempt to reproduce this in Mathematica is:

url = "http://www.massey.ac.nz/~pscowper/ts/cbe.dat";
CBE = Import[url];
elecTS = TemporalData[CBE[[2 ;; -1, 3]], {{1958}, Automatic, "Month"}];
DateListPlot[elecTS["Path"], Joined -> True, AspectRatio -> 0.2]

How do I continue this to perform the decomposition wholly within Mathematica?

Answer 1

Based on @b.gatessucks answer and on @RahulNarain comment tip, I created this functions for the multiplicative decompose case. I changed @b.gatessucks method for seasonality to keep it closer from R method, and used TemporalData to easily handle time interval.

decompose[data_,startDate_]:=Module[{dateRange,plot,plotOptions,observedData,observedPlot,trendData,trendPlot,dataDetrended,seasonalData,seasonalPlot,randomData,randomPlot},

    dateRange={{startDate,1},Automatic,"Month"};

    (*Setting Plot Options*)
    plotOptions={AspectRatio -> 0.2,ImageSize-> 600,ImagePadding->{{30,10},{10,5}},Joined->True};

    (*Observed data*)
    observedData=TemporalData[data,dateRange]["Path"];
    observedPlot=DateListPlot[observedData,Sequence@plotOptions];

    (*Extracting trend component*)
    trendData=Interpolation[MovingAverage[observedData,12],InterpolationOrder->1];
    trendPlot=DateListPlot[{#,trendData[#]}&/@observedData[[All,1]],Sequence@plotOptions]//Quiet;
    dataDetrended=N@{#[[1]],#[[2]]/trendData[#[[1]]]}&/@observedData//Quiet;

    (*Extracting seasonal component*)
    seasonalData=Mean/@Flatten[Partition[N@dataDetrended[[All,2]],12,12,1,{}],{2}];
    seasonalData=TemporalData[PadRight[seasonalData,Length[data],seasonalData],dateRange]["Path"];
    seasonalPlot=DateListPlot[seasonalData,Sequence@plotOptions];

    (*Extracting ramdon componet*)
    randomData=TemporalData[dataDetrended[[All,2]]/seasonalData[[All,2]],dateRange]["Path"];
    randomPlot=DateListPlot[randomData,Sequence@plotOptions];


    (*Ploting data*)
    plot=Labeled[
            Grid[Transpose[{Rotate[Style[#,15,Bold],90\[Degree]]&/@{"observed","trend","seasonal","random"}
                          ,{observedPlot,trendPlot,seasonalPlot,randomPlot}}
                          ]
                ]
            ,Style["Decomposition of multiplicative time series",Bold,17]
            ,Top
            ]

]

Using the functions like this:

rawData = Import["http://www.massey.ac.nz/~pscowper/ts/cbe.dat"][[2 ;;, 3]];
decompose[rawData, 1958]

We get:

Almost exactly as in R!

I say "almost" because R don't use interpolation, so the MoveAverage lost 12 point in R that we don't lose in this function due to interpolation method. I prefer to keek the ticks in each plot, I find it's better to read. It's a question of personal options.

Answer 2

While you come back with a version 9 solution here is an old school approach :

The first entry is labels so I removed it :

rawData = Import["http://www.massey.ac.nz/~pscowper/ts/cbe.dat"][[2 ;;]];

Added the dates to the imported data : they are monthly dates starting {1958, 1, 1} :

data = Transpose[{NestList[DatePlus[#, {1, "Month"}] &, {1958, 1, 1}, Length[rawData] - 1], 
                  rawData[[All, 3]]}];

plotData = DateListPlot[data, Joined -> True, PlotLabel -> "Data"];

The trend is a 12 - point moving average, found by trial and error and common sense :

mA = 12;
trend = Interpolation[Transpose[{AbsoluteTime[#] & /@ data[[1 ;; -mA, 1]], 
                      MovingAverage[data[[All, 2]], mA]}], InterpolationOrder -> 0];

plotTrend =  DateListPlot[{#, trend[AbsoluteTime[#]]} & /@ data[[1 ;; -mA, 1]], PlotLabel -> "Trend"];

Since we're doing a multiplicative decomposition we divide the starting data by its trend :

dataDetrended = {#[[1]], #[[2]]/trend[AbsoluteTime[#[[1]]]]} & /@ data[[1 ;; -mA]];

Seasonality is just a sine with a 1-year period :

seasonality = NonlinearModelFit[{AbsoluteTime[#[[1]]], #[[2]]} & /@ dataDetrended, 
                 a + b Sin[2 \[Pi] t/(365.25 86400 ) + f], {a, b, f}, t];

plotSeasonality = DateListPlot[{#, seasonality[AbsoluteTime[#]]} & /@ data[[1 ;; -mA, 1]], Joined -> True, PlotLabel -> "Seasonality"];

The random part is what is left after diving the de-trended data by the seasonal factor.

random = {#[[1]], #[[2]]/seasonality[AbsoluteTime[#[[1]]]]} & /@ dataDetrended[[1 ;; -mA]];

plotRandom = DateListPlot[random, Joined -> True, PlotLabel -> "Random"];

Final result :

GraphicsColumn[{plotData, plotTrend, plotSeasonality, plotRandom}, ImageSize -> 350]

Answer 3

This looks like very useful function for e.g. economic time series. If this can be done easily in R does it need to be done in Mathematica? With the qualifier that this is my first day testing 9 he is my attempt at the alternative:

Needs["RLink`"]
InstallR[];
data = REvaluate["{
      url <- \"http://www.massey.ac.nz/~pscowper/ts/cbe.dat\"
      CBE <- read.table(url, header = T)
      Elec.ts <- ts(CBE[, 3], start = 1958, freq = 12)
      decompose(Elec.ts, type = \"multi\")
   }"]

From here you can extract the decomposed data:

Normal@TemporalData[data[[1, x, 1]], {{1958}, Automatic, "Month"}]

where x=1, 2, 3, 4 give the observed, seasonal, trend, and random data. A solely Mathematica based approach would presumably require EstimatedProcess and SARIMAProcess as per some of the docs examples but I haven't got that far in learning 9.