Time-series decomposition in Mathematica
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 random component*)
randomData=TemporalData[dataDetrended[[All,2]]/seasonalData[[All,2]],dateRange]["Path"];
randomPlot=DateListPlot[randomData,Sequence@plotOptions];
(*Plotting 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 MovingAverage
lost 12 point in R that we don't lose in this function due to interpolation method.
I prefer to keep the ticks in each plot, I find it's better to read. It's a question of personal options.
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]
In my experience economists tend not to write time-series decomposition algorithms of the their own but rather use well established ones.
Among them two most well known are ARIMA-X13 and TRAMO-SEATS. Both are implemented by the US Census Bureau and executables are available here.
I've tried and implemented a simple package that calls the CB's files in the background and returns seasonally adjusted series. It very much follows the strategy that I learned by reading this answer on exporting multipage pdfs in mma on SE.
To use ARIMA-X13 in mma you should:
- Download Economica package here and put it into your $UserBaseDirectory.
- Download X13AS.EXE file from the CB's site and put it into "C:\Windows\System32"
Then you might extract the trend component of your series:
unemp = Rest@
QuandlData[
"ILOSTAT/UNE_TUNE_NB_SEX_T_AGE_AGGREGATE_TOTAL_Q_RUS", {{2001, 1,
1}, {2014, 1, 1}}]
unempSA = SeasonalAdjustmentX13[Sort@unemp]
DateListPlot[{unemp, unempSA}]
ARIMA-X13 allows also to extract cyclical component and noise, however I've yet to add them to that package.