The Season Package (by Ekkehart Schlicht)

Revision 2 of June 2006. Note that the output format of Season[ ] has changed!

The package implements Schlicht's (1984) seasonal adjustment method. It decomposes a time series \(TraditionalForm\`x\_t, \ \(\(t\)\(\ \)\(∈\)\(\ \)\({1,  ... , T}\)\(\ \)\)\)into a trend \(TraditionalForm\`y\_t, \ \(\(t\)\(\ \)\(∈\)\(\ \)\({1,  ... , T}\)\(\ \)\)\), a seasonal component \(TraditionalForm\`z\_t, \ \(\(t\)\(\ \)\(∈\)\(\ \)\({1,  ... , T}\)\(\ \)\)\), and an irregular component \(TraditionalForm\`u\_t, \ \(\(t\)\(\ \)\(∈\)\(\ \)\({1,  ... , T}\)\(\ \)\)\)such that (Subscript[x, t]) = 
        (Subscript[y, t]) + (Subscript[z, t]) + (Subscript[u, t]) for all  t. It combines the trend filter proposed by Leser (1961) (also known as the Hodrick-Prescott filter), the seasonal filter proposed by Schlicht and Pauly (1983) and the orthogonal parametrization proposed by Schlicht (1984). In contrast to prevailing methods, it is based on an explicit statistical model (state-space). 
Season is appropriate for smoothing time series with higher than annual frequency.  Season is inappropriate for smoothing annual data or deseasonalized time series; for this, use the HPFilter package by Johannes Ludsteck.   
This package provides the following functions for seasonal adjustment:

The method needs the time series x and the length of the season s as an input (e.g. s=4 for quarterly data or s=12 for monthly data). The function Season[x,s] computes smoothing parameters α and γ by a maximum likelihood method and splits the time series, using these weights. Alternatively, the parameters can be provided in Season[x,s,α,γ]. The parameter α controls the smoothness of the trend - the larger α, the smoother the trend. The parameter γ controls the rigidity of the seasonal pattern - the larger γ, the more time-invariant is the seasonal pattern. Both parameters interact: Increasing α while keeping γ constant will shift some variabilty from the trend to the seasonal component, for example (and the rest to the irregular component). 
The method assumes that trend, season and irregular component are additive components of the time series. For some time series, such as trended time series, a multiplicative formulation is preferable. In this case, take the logarithms of the original data as the input.

This loads the package

<< Season`Season`

Here is a series with 128 elements - the time series of monthly Gereman unemployment rates for 1992:01 through 2002:08.

x={8.3, 8.2, 7.7, 7.6, 7.4, 7.3, 7.8, 7.7, 7.5, 7.6, 7.7, 
            8.1, 9, 9, 8.7, 8.6, 8.4, 8.5, 9.1, 9.1, 8.9, 9.1, 9.2, 9.6, 
            10.5, 10.5, 10.1, 9.9, 9.5, 9.3, 9.6, 9.4, 9.1, 8.9, 8.9, 9.2, 
            10, 10, 9.6, 9.4, 9, 9, 9.4, 9.3, 9.2, 9.2, 9.3, 9.9, 10.8, 11.1, 
            10.8, 10.4, 10, 9.9, 10.2, 10.2, 10.1, 10.1, 10.3, 10.8, 12.2, 
            12.2, 11.7, 11.3, 11.1, 11, 11.4, 11.4, 11.2, 11.2, 11.3, 11.8, 
            12.6, 12.6, 12.1, 11.4, 10.9, 10.5, 10.7, 10.6, 10.3, 10.1, 10.2, 
            10.9, 11.5, 11.6, 11.1, 10.7, 10.2, 10.1, 10.3, 10.3, 10.1, 9.9, 
            10, 10.3, 11, 10.9, 10.6, 9.8, 9.3, 9.1, 9.3, 9.3, 9, 8.9, 8.9, 
            9.3, 10, 10.1, 9.8, 9.5, 9, 8.9, 9.2, 9.2, 9, 9, 9.2, 9.6, 10.4, 
            10.4, 10, 9.7, 9.5, 9.5, 9.7, 9.6};

The following code calls the decomposition routine for monthly data with smoothing parameter values α=100 and γ=100 and returns the list {{x,y,z,u},{α,γ}}.

{{x, y, z, u}, {a, g}, varU} = 
        Season[x, 12, 100,  100];

Here is a plot of the original time series x (red), the trend y (blue), and the seasonal component z (green).

<< Graphics`MultipleListPlot`
MultipleListPlot[x, y, z , PlotJoined -> True, 
        PlotStyle -> {{RGBColor[1, 0, 0]}, {RGBColor[0, 0, 
                1]}, {RGBColor[0, 1, 0]}}, 
        SymbolShape -> None];
Print["Smoothing parameters: α = ", 
        a, ",  γ = ", g];
\(\(Print[\*""Implied variances: \!\(σ\_u\^2\) = "", \ 
        varU, \ \*"", \!\(σ\_v\^2\) = "", 
        varU/a, \*"", \!\(σ\_w\^2\) = "", 
        varU/g];\)
\)

If Season is called without smoothing parameters, it computes the optimum smoothing parameters by a maximum-likelihood method and returns the list  {{x,y,z,u},{α,γ},(σ^2)}.

{{x, y, z, u}, {a, g}, var} = Season[x, 12];

Here is a plot of the original time series x (red), the trend y (blue), and the seasonal component z (green). It can be seen that the smoothing parameters are estimated smaller than the parameters used in the previous example. In particular, the parameter α for the trend is estimated as being quite small, entailing a rougher trend.

<< Graphics`MultipleListPlot`
MultipleListPlot[x, y, z , PlotJoined -> True, 
        PlotStyle -> {{RGBColor[1, 0, 0]}, {RGBColor[0, 0, 
                1]}, {RGBColor[0, 1, 0]}}, 
        SymbolShape -> None];
Print["Smoothing parameters: α = ", 
        a, ",  γ = ", g];
\(\(Print[\*""Estimated variances: \!\(σ\_u\^2\) = "", \ 
        varU, \ \*"", \!\(σ\_v\^2\) = "", 
        varU/a, \*"", \!\(σ\_w\^2\) = "", varU/g];\)\)

Here is an example:

{{x, y, z, u}, {a, g},  var} = 
        Season[x, 
          12,  "SearchRange"  ->  {{1, 2}, {1, 
                2}}];

Here is an example:

{{x, y, z, u}, {a, g},  var} = 
        Season[x, 12, AccuracyGoal -> 8, 
          MaxIterations -> 5];

Here is an example:

Likelihood[x, 12, 2, 3]

Here is an example:

LLPlot[x, 12, {0.9, 0.95}, {1.8, 2}, 
      PlotPoints  -> 20]

Evaluate the following code to clear symbols

Clear[ x, y, z, u, a, g];