$\mathrm{path}:=\mathrm{FileTools}:-\mathrm{JoinPath}\left(\left[\mathrm{kernelopts}\left(\'\mathrm{datadir}\'\right)\,''datasets''\,''air\_passengers.csv''\right]\right)\:$

$\mathrm{data}:=\mathrm{ImportMatrix}\left(\mathrm{path}\right)$

The data has one column of dates (year) and one column of data (monthly passengers). The first row is a header.

$\mathrm{data}\left[1..5\right]$

To work with this data, construct a TimeSeries object. Such an object can contain one or more data sets, measured at a common set of time points, as well as data headers and other metadata. You can construct this time series object as follows:

$\mathrm{ts}:=\mathrm{TimeSeries}\left(\mathrm{data}\,\'\mathrm{header}\'\=1\,\'\mathrm{dates}\'\=1\,\'\mathrm{period}\'\=12\right)$

The options above indicate that the first row and column contain a header and the dates, respectively, and that you expect any seasonal characteristics to occur with period 12.

To inspect the data, you can use the GetData, GetDates, and GetHeaders commands.

$\mathrm{GetData}\left(\mathrm{ts}\right)$

$\mathrm{GetDates}\left(\mathrm{ts}\right)$

$\mathrm{GetHeaders}\left(\mathrm{ts}\right)$

$\left[''Monthly\; Passengers''\right]$

Alternatively, you can plot the data using the TimeSeriesPlot command.

$\mathrm{TimeSeriesPlot}\left(\mathrm{ts}\right)$

We can also look for seasonal trends in our data using the SeasonalSubseriesPlot command. The following plot shows the number of passengers on a monthly basis.

$\mathrm{SeasonalSubseriesPlot}\left(\mathrm{ts}\,\mathit{seasonnames}\=\left[''January''\,''February''\,''March''\,''April''\,''May''\,''June''\,''July''\,''August''\,''September''\,''October''\,''November''\,''December''\right]\,\mathrm{size}\=\left[800\,300\right]\right)$

You can now have Maple select a suitable model from a family of 30 related models, and adjust it to this time series. In fact, you will adjust it to the first 10 years of data; this can be done by specifying the time range as an index.

$\mathrm{tstrain}≔\mathrm{ts}\left[..''1958-12''\right]$

$\mathrm{tsverify}≔\mathrm{ts}\left[''1959'' ..\right]$

The family of models used is an exponential smoothing model.

$\mathrm{model}:=\mathrm{ExponentialSmoothingModel}\left(\mathrm{tstrain}\right)$

$\mathrm{model}:=\mathrm{<\; an\; ETS(M,M,M)\; model\; >}$

You can observe that the best model in this case has multiplicative errors, a multiplicative trend, and multiplicative seasonal information. You can get the details on this exponential smoothing model using the GetParameters command.

$\mathrm{GetParameters}\left(\mathrm{model}\right)$

Use the model to predict two years of future data using the forecast command: the actual data for this time range is in $\mathrm{tsverify}$.

$\mathrm{forecast}:=\mathrm{Forecast}\left(\mathrm{model}\&comma;\mathrm{tstrain}\&comma;24\right)$

The forecast is itself a time series object that can be inspected using GetData, GetHeaders, and other commands.

$\mathrm{GetHeaders}\left(\mathrm{forecast}\right)$

$\left[''Monthly\; Passengers\; (forecast)''\right]$

$\mathrm{TimeSeriesPlot}\left(\mathrm{tstrain}\&comma;\mathrm{tsverify}\&comma;\mathrm{forecast}\right)$

There is reasonable agreement between the forecast and the verification data. You get a better handle on this if you include confidence intervals for all data points. You can then see how often the true data falls within the boundaries of the confidence intervals.

$\mathrm{confidence}:=\mathrm{Forecast}\left(\mathrm{model}\&comma;\mathrm{tstrain}\&comma;24\&comma;\mathrm{output}\&equals;\mathrm{confidenceintervals}\left(80\&comma;95\right)\right)$

$\mathrm{TimeSeriesPlot}\left(\mathrm{tstrain}\&comma;\mathrm{forecast}\&comma;\mathrm{confidence}\&comma;\left[\mathrm{tsverify}\&comma;\mathrm{color}\&equals;''Red''\&comma;\mathrm{thickness}\&equals;3\right]\right)$

Use the model to decompose the data set into several components; the number of components depends on the model. This is done using the Decomposition command. For an exponential smoothing model, there are always level and residual components. There may also be trend and seasonal components depending on whether or not the model has those properties. This can be used to, for example, correct for seasonal influences or smooth data.

$\mathrm{decomposition}:=\mathrm{Decomposition}\left(\mathrm{model}\&comma;\mathrm{tstrain}\right)$

$\mathrm{TimeSeriesPlot}\left(\mathrm{tstrain}\&comma;\mathrm{decomposition}\&comma;\mathrm{split}\&equals;\mathrm{pertimeseries}\right)$

Because the seasonal component is multiplicative, you can compensate for it by dividing the original data by it.

$\mathrm{trainingdata}:=\mathrm{GetData}\left(\mathrm{tstrain}\right)$

$\mathrm{seasonaldata}:=\mathrm{GetData}\left(\mathrm{decomposition}\left[..\&comma; ''Monthly\; Passengers\; (seasonal)''\right]\right)$

$\mathrm{nonseasonaldata} ≔ \mathrm{TimeSeries}\left(\mathrm{trainingdata} /~ \mathrm{seasonaldata}\&comma; \&apos;\mathrm{header}\&apos; \&equals; ''Monthly\; Passengers\; (deseasonalized)''\&comma; \&apos;\mathrm{dates}\&apos; \&equals; \mathrm{GetDates}\left(\mathrm{tstrain}\right)\right)$

$\mathrm{TimeSeriesPlot}\left(\mathrm{nonseasonaldata}\right)$

$\mathrm{merged}:=\mathrm{Join}\left(\mathrm{tstrain}\&comma;\mathrm{forecast}\right)$

$\mathrm{TimeSeriesPlot}\left(\mathrm{merged}\right)$

Another command, Difference, applies a differencing transformation in a flexible way. The LogTransform command applies a logarithm transformation, and BoxCoxTransform generalizes that to a general Box-Cox transformation.

$\mathrm{differenced}:=\mathrm{Apply}\left(\mathrm{Difference}\&comma;\mathrm{ts}\right)$

$\mathrm{logs}:=\mathrm{Apply}\left(\mathrm{LogTransform}\&comma;\mathrm{ts}\right)$

$\mathrm{boxcox}:=\mathrm{Apply}\left(\mathrm{BoxCoxTransform}\left(\mathrm{\&lambda;}\&equals;\frac{1}{3}\right)\&comma;\mathrm{ts}\right)$

$\mathrm{TimeSeriesPlot}\left(\mathrm{ts}\&comma;\mathrm{differenced}\&comma;\mathrm{logs}\&comma;\mathrm{boxcox}\&comma;\mathrm{split}\&equals;\mathrm{pertimeseries}\right)$

$\mathrm{general\_model}:=\mathrm{ExponentialSmoothingModel}\left(\&apos;\mathrm{seasonal}\&equals;\left\{A\&comma;M\right\}\&apos;\right)$

$\mathrm{general\_model}:=\mathrm{<\; an\; ETS(*,*,*)\; model\; >}$

You can now specialize this to all individual model formulations represented by the general model. That is, where the seasonal component is additive or multiplicative. Some models are excluded by default because they are subject to numerical difficulties: the forecasts have infinite variance. (They can be included by overriding an option to Specialize.)

$\mathrm{individual\_models}:=\mathrm{Specialize}\left(\mathrm{general\_model}\&comma;\mathrm{ts}\right)$

$\mathrm{individual\_models}:=\left[\mathrm{<\; an\; ETS(A,A,A)\; model\; >}\&comma;\mathrm{<\; an\; ETS(A,Ad,A)\; model\; >}\&comma;\mathrm{<\; an\; ETS(A,N,A)\; model\; >}\&comma;\mathrm{<\; an\; ETS(M,A,A)\; model\; >}\&comma;\mathrm{<\; an\; ETS(M,A,M)\; model\; >}\&comma;\mathrm{<\; an\; ETS(M,Ad,A)\; model\; >}\&comma;\mathrm{<\; an\; ETS(M,Ad,M)\; model\; >}\&comma;\mathrm{<\; an\; ETS(M,M,M)\; model\; >}\&comma;\mathrm{<\; an\; ETS(M,Md,M)\; model\; >}\&comma;\mathrm{<\; an\; ETS(M,N,A)\; model\; >}\&comma;\mathrm{<\; an\; ETS(M,N,M)\; model\; >}\right]$

Each individual model now has a slightly different set of model equations. They all still have a number of parameters and initial values for the model's state variables and need to be optimized for the best fit. This optimization process consists of a number of iterations. In each iteration, Maple picks a new set of values for the parameters and initial state values, then runs the simulation using the model, and finally, computes the deviations from the actual observed data. The first set of values (which initializes the optimization process) is computed directly from the actual data, by the Initialize command.

$\mathrm{initialization\_tables}:=\mathrm{map}\left(\mathrm{Initialize}\&comma;\mathrm{individual\_models}\&comma;\mathrm{ts}\right)\&colon;$

Here is an example of these initialization values:

${\mathrm{individual\_models}}_1\&comma;{\mathrm{initialization\_tables}}_1$

You can now perform the optimization.

$\mathrm{map}\left(\mathrm{Optimize}\&comma;\mathrm{individual\_models}\&comma;\mathrm{ts}\right)$

$\left[-579.1189017\&comma;-587.8357482\&comma;-587.5947092\&comma;-577.5718572\&comma;-527.4987042\&comma;-559.4653237\&comma;-528.6819747\&comma;-527.8440927\&comma;-527.3996887\&comma;-565.0019522\&comma;-535.5928547\right]$

The optimization process sets all parameters and initial state values to the optimal values found. The Optimize command returns a measure for how close the fit between the simulation run and the actual data is. However, these models all have different numbers of parameters. If a similar fit is achieved with a model that has fewer parameters, then the principle of parsimony says you should prefer the latter model. This is quantified in a so-called information criterion: a function that takes both the closeness of the fit and the number of parameters into account. The TimeSeriesAnalysis package has three information criteria built in: two versions of Akaike's Information Criterion (AICc, which includes a correction for small data sizes, and AIC, which does not), and the Bayesian Information Criterion (BIC). Let us compare the BIC for all these models.

$$\begin{gathered} \mathbf{for}m\mathbf{in}\mathrm{individual\_models}\mathbf{do} \\ \mathrm{print}\left(m\&comma;\mathrm{BIC}\left(m\&comma;\mathrm{ts}\right)\right) \\ \mathbf{end} \mathbf{do}\&colon; \end{gathered}$$

$\mathrm{<\; an\; ETS(A,A,A)\; model\; >}\&comma;1412.79583526645$

$\mathrm{<\; an\; ETS(A,Ad,A)\; model\; >}\&comma;1445.84094061076$

$\mathrm{<\; an\; ETS(A,N,A)\; model\; >}\&comma;1422.76337107009$

$\mathrm{<\; an\; ETS(M,A,A)\; model\; >}\&comma;1722.03500867297$

$\mathrm{<\; an\; ETS(M,A,M)\; model\; >}\&comma;1133.89411113851$

$\mathrm{<\; an\; ETS(M,Ad,A)\; model\; >}\&comma;1945.38209291183$

$\mathrm{<\; an\; ETS(M,Ad,M)\; model\; >}\&comma;1141.84855312639$

$\mathrm{<\; an\; ETS(M,M,M)\; model\; >}\&comma;1139.65948604423$

$\mathrm{<\; an\; ETS(M,Md,M)\; model\; >}\&comma;1141.03566969852$

$\mathrm{<\; an\; ETS(M,N,A)\; model\; >}\&comma;2087.66988840404$

$\mathrm{<\; an\; ETS(M,N,M)\; model\; >}\&comma;1150.95261586989$

From the results, the $\left(M\&comma;A\&comma;M\right)$, $\left(M\&comma;\mathrm{Ad}\&comma;M\right)$,$\left(M\&comma;M\&comma;M\right)$, $\&lpar;M\&comma;N\&comma;M$), and ($M\&comma;\mathrm{Md}\&comma;M\&rpar;$models give the best results in terms of the Bayesian information criterion.

Box, G.E.P., Jenkins, G.M., and Reinsel, G.C. (1976) Time Series Analysis, Forecasting and Control. Third Edition. Holden-Day. Series G.

Hyndman, R.J. and Athanasopoulos, G. (2013) Forecasting: principles and practice. http://otexts.org/fpp/. Accessed on 2013-10-09.

Hyndman, R.J., Koehler, A.B., Ord, J.K., and Snyder, R.D. (2008) Forecasting with Exponential Smoothing: The State Space Approach. Springer Series in Statistics. Springer-Verlag Berlin Heidelberg.