An analysis of the Lynx data set

We model the celebrated data recorded for the number of lynx trapped annually in northwest Canada from 1821 to 1934. We first take the Log of the data, lynxdata, and remove its mean.  The result is called ZeroMeanLoglynxdata.

PlotZeroMeanLoglynxdata

Here are the correlation function and partial correlation function for the Lynx data set out to a lag of 30.

DisplayTogetherArray[
      ListPlot[CorrelationFunction[ZeroMeanLoglynxdata,  30], options1], 
      ListPlot[PartialCorrelationFunction[ZeroMeanLoglynxdata,  30], 
        options2]];

The Hannan-Rissanen method is used to select from among 30 candidate models.

models = HannanRissanenEstimate[ZeroMeanLoglynxdata, 15, 15, 15, 30];

We compute the AIC's of models, choose the one with the smallest value, and assign it to the variable initialModel.

Determine intialModel

initialModel is an ARMA(11,1) model. We will attempt to refine initialModel by using it as input to the conditional maximum likelihood estimate. We assign this model to the variable refinedModel.

refinedModel = 
        ConditionalMLEstimate[ZeroMeanLoglynxdata, 
          initialModel];

Indeed, refinedModel has a smaller PortmanteauStatistic than initialModel:

PortmanteauStatistic[Residual[ZeroMeanLoglynxdata, #], 
          30] & /@ {initialModel, refinedModel}

The correlation and partial correlation functions of refinedModel are quite similar to those of the original data.

Display Correlation and Partial Correlation Functions

We predict the next 15 values using BestLinearPredictor. The predicted data are show in red.

{xhat, error} = BestLinearPredictor[ZeroMeanLoglynxdata, refinedModel, 15];
PlotData[Loglynxdata,xhat+μ];