(* ::Package:: *)

(* :Title: Season *)

(* :Author: Ekkehart Schlicht *)

(* :Summary:
   This package provides functions to decompose a given time series
   into trend, sesonal component and irregular component as described 
   Ekkehart Schlicht (1984). 
*)

(* :Mathematica Version: 6.0.3.0
	Runs also on Mathematica 8 to 10
	Program version 1.02 of October 10, 2005
	revised July 10, 2006, February 11, 2017, December 2020 
    

	
*)

BeginPackage["Season`Season`"];

Unprotect[{Season, LL, LLPlot}];


Season::usage = "Season[x,s,options] decomposes the time series x into trend, \
\
seasonal seasonal component, and irregular component according to the method \
\
by Schlicht (1984). The integer s denotes the length of the seasonal pattern \
\
(e.g. s=4 for quarterly data), s=12. for monthly data, or s=7 for daily data with \
\
a weekly pattern.

In case of success the function returns a list {{x,y,z,u},{alpha,gamma},varU} \
\
where x is the original time series, y is the trend, z is the seasonal component, \
\
and u is the irregular component. The smoothing constants used are alpha (for the \
trend) and gamma (for the seasonal component). varU gives the estimated variance \
of the irregular component. The variance of the trend disturbance is varU/alpha, \
the variance of the seasonal disturbance is varU/gamma.

If Season[x,s,alpha,gamma,options] is called, the smoothing constants are \
taken from the input rather than being computed.

Possible options for Season are all NMaximize options (except the Method \
option).

Furthermore it is possible to specify limits of the search range for the \
\
smoothing constants alpha and gamma using the option SearchRange -> {{aMin, aMax},{gMin,gMax}}. \
\
This option may also be used to fix one smoothing constant at some value \
\
and have the other one estimated. If a boundary binds, a warning is issued."

LL::usage = "LL[x,s,alpha,gamma] returns the log likelihood \
of \
{x,s,alpha,gamma}."

LLPlot::usage = "LLPlot[x,s,{aMin,aMax},{gMin,gMax}] plots the log likelihood \
\
in the range [aMin,aMax]x[gMin,gMax]. If options are given by the user, they \
are
handed over to Plot3D which is called by LLPlot."
 

Options[Season] = {SearchRange -> {{0.1, 1000}, {0.1, 1000}},
                   Method -> {"NelderMead", "PostProcess" -> False}};
 
Begin["`Private`"];

makeP[T_] := SparseArray[
    {{i_, i_} -> 1,{i_, j_} /; j == i + 1 -> -2,{i_, j_} /; 
    j == i + 2 -> 1},{T - 2, T}];
    
makeR[T_,s_] := SparseArray[
    {{i_, j_} /; (j <=(i + s-1))&&(j>=i) -> 1},
    {T - s+1, T}];

makeQ[T_,s_] := (1/(s-1)) SparseArray[
    {{i_, j_} /; (j < (i + s-1))&&(j>=i) -> (j-i+1)},
    {T - s+1, T-1}]; 
    
init[x_,s_]:=
  Block[{T,im,p,q,r,pp,ppI,qq,qqI,rri,rqqIr,Sy,Sz},
    T=Length[x]; im = N[IdentityMatrix[T]];
    p = N[makeP[T]]; q = N[makeQ[T,s]]; r = N[makeR[T,s]];
    pp = Transpose[p].p; ppI = Inverse[p.Transpose[p]];
    qq = q.Transpose[q]; qqI = Inverse[qq];
    rrI = Inverse[r.Transpose[r]];
    rqqIr = Transpose[r].qqI.r;
    Sy = Transpose[p].ppI.ppI.p;
    Sz = Transpose[r].rrI.qq.rrI.r;
    {im,Sy,Sz,pp,rqqIr}];
   
split[x_,s_,{im_,Sy_,Sz_,pp_,rqqIr_},a_?NumberQ,g_?NumberQ] := Module[
  {ldet, H, H1, H2, sol, y, sq, crit, T = Length[x]},
    ldet=Log[Det[im + Sy/a + Sz/g]];      
    H=SparseArray[ArrayFlatten[{{a pp + im, im},{im, g rqqIr + im}}]];
    H=SparseArray[(1/2) (H + Transpose[H])];
    sol=LinearSolve[H, Join[x,x], Method->Cholesky];
    y = Take[sol,T]; z = Drop[sol,T]; u = x-y-z;
    sq = u.u + a y.pp.y + g z.rqqIr.z;
    {{x,y,z,u}, {a,g},sq/T, -ldet/T - Log[sq]}];

    
setOptions[opts___Rule]:= Block[
  {userOpts, restOpts, allOpts, aMin, aMax, gMin, gMax},
  
  userOpts = DeleteCases[{opts},Rule[Method,_]];  
  restOpts = Complement[Options[Season], Sequence @@ List /@ userOpts];
  allOpts = Join[userOpts, restOpts];
  {{aMin, aMax}, {gMin,gMax}} = SearchRange /. allOpts;
  {allOpts, {{aMin, aMax}, {gMin, gMax}}}];

  
checkArgs[x_, s_]:= Which[
  s <= 1, Message[Season::"freqmin"]; Abort[],
  s > Length[x], Message[Season::"freqmax"]; Abort[]];

LH[x_,s_,mat_, a_?NumberQ,g_?NumberQ]:= Last[split[x,s,mat,a,g]];

Season[x_?VectorQ,s_?IntegerQ,opts___Rule]:=
  Module[{a,g, matrices},  

   checkArgs[x,s];
   matrices = init[x,s]; 
   {allOpts, {{aMin, aMax}, {gMin, gMax}}} = setOptions[opts]; 
   
    maxLH=Check[NMaximize[{LH[x,s,matrices,a,g],
           aMin <= a <= aMax, gMin <= g <= gMax}, {a,g}, 
           Sequence@@FilterRules[{allOpts}, Options@NMaximize]],
           $Failed];

   If[maxLH === $Failed,
      Message[Season::"maxfail"]; Abort[]]; (* else *)
           
    {a, g} = {a, g}/.Last[maxLH];
    
    If[((aMax-a)<0.001 || (a-aMin)<0.001), Message[Season::"cornera", a]];
    If[((gMax-g)<0.001 || (g-gMin)<0.001), Message[Season::"cornerg", g]];
    
    
    Most[split[x,s,matrices,a,g]]];
    
    
Season[x_?VectorQ,s_?IntegerQ,a_,g_]:=
  Block[{matrices},
    checkArgs[x,s];
    matrices = init[x,s];
    Most[split[x,s,matrices,a,g]]];
    
LL[x_?VectorQ,s_?IntegerQ,a_,g_] := (
	checkArgs[x,s];
   	matrices=init[x,s]; T=Length[x];   	
 	(LH[x,s,matrices,a,g] T + T Log[T] + T)); 
   
LLPlot[x_?VectorQ,s_?IntegerQ,{aMin_,aMax_},{gMin_,gMax_},opts___]:=
  	Block[{a,g},
 	matrices=init[x,s]; T=Length[x];
  	Plot3D[(LH[x,s,matrices,a,g] T + T Log[T] + \
T),{a,aMin,aMax},{g,gMin,gMax},opts]];
  
    
Season::"spar" = "Both smoothing constants must be larger than zero - \
calculation aborted.";

Season::"freqmin" = "The minimal length of the seasonal pattern is two - \
calculation aborted.";

Season::"freqmax" = "The length of the seasonal pattern must not exceed " <>
                    "the length of the time series - calculation aborted.";
    
Season::"cornera" = "Estimate `1` of alpha appears to be a corner solution. Try to widen the search range.";
Season::"cornerg" = "Estimate `1` of gamma appears to be a corner solution. Try to widen the search range.";
   
Season::"maxfail" = "Aborting due to problems in numerical maximization \
procedure." <>
         "\nConsider fine-tuning NMaximize options. Try the option MaxIterations\[Rule]500, for example.";
    
Protect[{Season, LL, LLPLot}];

End[];

EndPackage[];







