Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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:

enter image description here

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]

enter image description here

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

share|improve this question
    
I would also really like to how to do this solely in Mathematica. My guess is that some of the Version 9 built in time-series functions will make it relatively simple. –  Cameron Murray Dec 21 '12 at 3:38
1  
While I have no doubt that Mathematica is meta tool and all that, the fact remains that specialized software do the job better-like for decomposition and econometric analysis, nothing beats Eviews –  rselva Dec 21 '12 at 8:31
    
@rselva Asking in the most constructive spirit, could you provide some examples illustrating your statement ? –  b.gatessucks Dec 21 '12 at 10:54
    
@rselva Agree on Eviews (I wish I could afford a personal license). FWIW, as much as I love MMA, using it for econometrics is at best a challenge, and at worse a time killer. I wish MMA's design could be made more flexible for econometrics, but I understand why this would be a messy chore-and-a-half. I plan to learn some R myself (and am pleased they've built in this new capability). For those interested in a version of R that is automatically linked up to Excel, check out: rcom.univie.ac.at –  telefunkenvf14 Dec 21 '12 at 20:57
    
This question would be easier to answer if you could find a reference to the algorithm that R's time series decomposition is using. It might even be described in the documentation. Asking to reverse-engineer the functionality of a magical third-party black box is still a valid question, but a much harder one (as evidenced by your comment on b.gatessucks's answer). –  Rahul Narain Dec 23 '12 at 1:30
show 3 more comments

3 Answers

up vote 12 down vote accepted

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:

enter image description here

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.

share|improve this answer
add comment

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]

enter image description here

share|improve this answer
    
Very cool! My only question is about the seasonal part. I believe that the Sin assumption is too specific, but I don't know yet a more general way to do that, maybe with more terms, like a fourrier approach. –  Murta Dec 22 '12 at 23:09
add comment

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.

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.