Canadian lynx series

A first time series model of the Canadian lynx data was fitted by P.A.P. Moran [a13] in 1953. He observed that the cycle is very asymmetrical with a sharp and large peak and a relatively smooth and small trough. The log transformation gives a series which appears to vary symmetrically about the mean. As the actual population of lynx is not exactly proportional to the number caught, a better representation would perhaps be obtained by incorporating an additional "error of observation" in the model, thereby resulting in a more complicated model. The log transformation substantially reduces the effect of ignoring this error of observation; therefore, after Moran, nearly all the time series analysis of the lynx data in the literature have used the log-transformed data. Let

$$ X _ {t} = $$

$$ = { \mathop{\rm log} } _ {10 } ( \textrm{ number recorded as trapped in year 1820 } + t ) $$

( $ t = 1 \dots 114 $). Because of the apparently slow damping in the amplitude of the sample correlogram, Moran discarded the idea of a sinusoidal-shape model and proposed an $ { \mathop{\rm AR} } ( 2 ) $- model.

In 1977, M.J. Campbell and A.M. Walker [a2] believed that an appropriate model of the lynx data should be, in some sense, "between" a pure harmonic model and a pure auto-regression. Subsequently this led them to combining a harmonic component with an $ { \mathop{\rm AR} } ( 2 ) $- model. Two models with frequencies $ 9.5 $, $ { \mathop{\rm CW} } ( 9.5 ) $, and $ 9.63 $, $ { \mathop{\rm CW} } ( 9.63 ) $, were recommended. At about the same time, an $ { \mathop{\rm AR} } ( 11 ) $- model based on the Akaike information criterion was fitted [a20]. In the discussion of the above two papers, D.R. Cox [a4] suggested a polynomial model.

In 1979, R.J. Bhansali [a1] used a mixed spectrum analysis to analyze the lynx data.

Using the Canadian lynx data as a case study, H. Tong and K.S. Lim [a24] fitted a class of non-linear models called the self-exciting threshold auto-regressive model (SETAR model) to the log-transformed lynx data. They demonstrated that this model has interesting features in non-linear oscillations, such as jump resonance, amplitude-frequency dependency, limit cycles, subharmonics, and higher harmonics. Later, in 1981, it was discovered that the self-exciting threshold auto-regressive model also generates chaos [a11]. In the discussion of [a24], T. Subba Rao and M.M. Gabr [a16] proposed a subset bilinear model, $ { \mathop{\rm SBL} } ( 11 ) $, to the first one hundred log-transformed lynx data and an $ { \mathop{\rm SBL} } ( 9 ) $ to the first one hundred original lynx data.

In 1981, they used the maximum-likelihood estimation coupled with the Akaike information criterion to fit a subset bilinear model $ { \mathop{\rm SBL} } ( 12 ) $ to the first one hundred log-transformed lynx data [a5]. Their model was able to produce small values of the noise variance and the mean-squared errors of the one-step-ahead predictions, but it failed to detect the inherited behaviour of the data.

V. Haggan and T. Ozaki [a8] fitted an exponential auto-regressive model, $ { \mathop{\rm EXPAR} } ( 11 ) $, another class of non-linear models, to the mean-deleted log-transformed lynx data in 1981. Ozaki [a15] felt that the almost symmetric series generated by this model was unsatisfactory. Subsequently, in 1982, he fitted two more exponential auto-regressive models to the full set of log-transformed lynx data with mean deleted. One of them, $ { \mathop{\rm EXPAR} } ( 2 ) $, could reproduce the asymmetric limit cycle structure of the lynx data; the other, $ { \mathop{\rm EXPAR} } ( 9 ) $, with smaller variance of fitted residuals, was believed to be more appropriate for forecasting.

D.F. Nicholls and D.G. Quinn [a14], in 1982, fitted another new class of time series models, called a random coefficient auto-regressive model (RCA model) to the first one hundred log-transformed lynx data, using a maximum-likelihood method or the conditional least-squares method.

In 1984, Haggan, S.M. Heravi and M.B. Priestley [a7] fitted a state-dependent model (SDM) of $ { \mathop{\rm AR} } ( 2 ) $ to the log-transformed lynx data.

Using the revised computer program in [a21], a $ { \mathop{\rm SETAR} } ( 2;5,2 ) $ was fitted to the first one hundred log-transformed lynx data, and a $ { \mathop{\rm SETAR} } ( 2;7,2 ) $ was fitted to the full set of log-transformed data as follows:

$$ X _ {t} = 0.546 + 1.032X _ {t - 1 } - 0.173X _ {t - 2 } + $$

$$ + 0.171X _ {t - 3 } - 0.431X _ {t - 4 } + 0.332X _ {t - 5 } - $$

$$ - 0.284X _ {t - 6 } + 0.210X _ {t - 7 } + \epsilon _ {t} ^ {( 1 ) } $$

if $ X _ {t - 2 } \leq 3.116 $, and

$$ X _ {t} = 2.632 + 1.492X _ {t - 1 } - 1.324X _ {t - 2 } + \epsilon _ {t} ^ {( 2 ) } $$

if $ X _ {t - 2 } > 3.116 $. Here, $ { \mathop{\rm var} } ( \epsilon _ {t} ^ {( 1 ) } ) = 0.0258 $, $ { \mathop{\rm var} } ( \epsilon _ {t} ^ {( 2 ) } ) = 0.0505 $( pooled variance equals $ 0.0360 $).

This model was able to describe the biological features of the Canadian lynx data such as:

1) its cyclical behaviour of about 9–10 years per cycle;

2) the rise periods exceed the descent periods in the cycles;

3) the delay parameter of $ 2 $ in $ X _ {t - 2 } \leq 3.116 $ is associated with the biological cycle that a Canadian lynx is fully grown in the autumn of its second year and births of kittens (1–4 per litter) take place about 63 days after breeding in March–April;

4) the threshold estimate, $ 3.116 $, lies in the vicinity of the anti-mode of the histogram of the lynx data, which implies that there is insufficient information in the data to model more precisely the functional form of the dynamics over the state space near the sample mean.

A comparative study of some of the above models was carried out in [a12]. The $ { \mathop{\rm SETAR} } $ models were ranked to be the best among the models considered. For an extensive discussion of the lynx data, see [a22].

B.Y. Thanoon (1988) [a18] fitted the two subset $ { \mathop{\rm SETAR} } $( $ { \mathop{\rm SSETAR} } $) models which produced limit cycles with two subcycles giving an average period of $ 9.5 $ years. He commented that the $ { \mathop{\rm SSETAR} } $ detected the inherited behaviour of the data better than the full $ { \mathop{\rm SETAR} } $ model in terms of the auto-covariance function.

In 1989, R.S. Tsay [a25] fitted a two-thresholds $ { \mathop{\rm SETAR} } ( 3;1,7,2 ) $ when proposing a new procedure for testing and building $ { \mathop{\rm TAR} } $ models. Around the same time, techniques from multivariate analysis were applied, [a23]; namely, the principal coordinate analysis and dendograms to twelve time series models reported in the literature.

In 1991, G.H. Yu and Y.C. Lin [a29] suggested a subset auto-regressive model, $ { \mathop{\rm SAR} } ( 1,3,9,12 ) $, to the log-transformed lynx data when proposing a method for selecting a best $ { \mathop{\rm SAR} } $ model automatically.

Applying the cross-validatory approach, in 1992, B. Cheng and Tong [a3] found the embedding dimension of the lynx data to be $ 3 $.

Using the lynx data as an example, in 1993, J. Geweke and N. Terui [a6] proposed a Bayesian approach for deriving the exact posterior distributions of the delay and threshold parameters.

In 1994, T. Teräsvirta [a17] fitted a logistic smooth transition auto-regressive model ( $ { \mathop{\rm LSTAR} } ( 11 ) $), which had a limit cycle of 77 years with eight subcycles of lengths 9 and 10 years. The lynx data has been used for the non-parametric identification of non-linear time series in selecting significant lags [a19]. The lynx data for 1821–1924 has been used [a27] to estimate $ f _ {m} ( \cdot ) $ and $ \lambda _ {m} ( \cdot ) $, the $ m $- step Lyapunov-like index, where $ \lambda _ {m} ( x ) = { {df _ {m} } / {dx } } $, and the last ten data to check the predicted values. The data was also used [a28] for subset selection in non-parametric stochastic regression. The subset of lags $ 1 $, $ 3 $, and $ 6 $ are selected from the original lynx data.

In 1995, C. Kooperberg, C.J. Stone and Y.K. Truong [a9] used their automatic procedure to estimate the mixed spectral distribution of the log-transformed lynx data and found the lynx cycle to be $ 9.5 $ years.

In 1996, D. Lai [a10] used a BDS statistic to test the residuals from the models of Moran, $ { \mathop{\rm SETAR} } ( 2;2,2 ) $, $ { \mathop{\rm SETAR} } ( 2;7,2 ) $, $ { \mathop{\rm SETAR} } ( 3;1,7,2 ) $, $ { \mathop{\rm EXPAR} } ( 11 ) $, $ { \mathop{\rm EXPAR} } ( 2 ) $, and Cox's polynomial model. He concluded that Tong's $ { \mathop{\rm SETAR} } ( 2;7,2 ) $ was found to be the best. C. Wong and R. Kohn [a26] used a Bayesian approach for estimating non-parametrically an additive auto-regressive model for the lynx data.

References

[a1]	R.J. Bhansali, "A mixed spectrum analysis of the Lynx data" J. Roy. Statist. Soc. A , 142 (1979) pp. 199–209 MR0547237
[a2]	M.J. Campbell, A.M. Walker, "A survey of statistical work on the McKenzie River series of annual Canadian lynx trappings for the years 1821–1934, and a new analysis" J. Roy. Statist. Soc. A , 140 (1977) pp. 411–431; discussion: 448–468
[a3]	B. Cheng, H. Tong, "On consistent non-parametric order determination and chaos" J. Roy. Statist. Soc. B , 54 (1992) pp. 427–449
[a4]	D.R. Cox, "Discussion of papers by Campbell and Walker, Tong and Morris" J. Roy. Statist. Soc. A , 140 (1977) pp. 453–454
[a5]	M.M. Gabr, T. Subba Rao, "The estimation and prediction of subset bilinear time series models with applications" J. Time Ser. Anal. , 2 (1981) pp. 153–171 MR0640211
[a6]	J. Geweke, N. Terui, "Bayesian threshold autoregressive models for nonlinear time series" J. Time Ser. Anal. , 14 (1993) pp. 441–454 MR1243574 Zbl 0779.62073
[a7]	V. Haggan, S.M. Heravi, M.B. Priestley, "A study of the application of state-dependent models in non-linear time series analysis" J. Time Ser. Anal. , 5 (1984) pp. 69–102 MR758580 Zbl 0555.62071
[a8]	V. Haggan, T. Ozaki, "Modelling non-linear random vibrations using an amplitude-dependent autoregressive time series model" Biometrika , 68 (1981) pp. 189–196 MR614955
[a9]	C. Kooperberg, C.J. Stone, Y.K. Truong, "Logspline estimation of a possibly mixed spectral distribution" J. Time Ser. Anal. , 16 (1995) pp. 359–388 MR1342682 Zbl 0832.62083
[a10]	D. Lai, "Comparison study of AR models on the Canadian lynx data: A close look at BDS statistic" Comm. Stat. Data Anal. , 22 (1996) pp. 409–423 MR1411579 Zbl 0875.62412
[a11]	K.S. Lim, "On threshold time series modelling" , Univ. Manchester (1981) (Doctoral Thesis (unpublished))
[a12]	K.S. Lim, "A comparative study of various univariate time series models for Canadian lynx data" J. Time Ser. Anal. , 8 (1987) pp. 161–176 Zbl 0608.62116
[a13]	P.A.P. Moran, "The statistical analysis of the Canadian lynx cycle. I: structure and prediction" Aust. J. Zool. , 1 (1953) pp. 163–173
[a14]	D.F. Nicholls, B.G. Quinn, "Random coefficient autoregressive models: an introduction" , Lecture Notes in Statistics , 11 , Springer (1982) MR0671255 Zbl 0497.62081
[a15]	T. Ozaki, "The statistical analysis of perturbed limit cycle processes using nonlinear time series models" J. Time Ser. Anal. , 3 (1982) pp. 29–41 MR0660394 Zbl 0499.62079
[a16]	T. Subba Rao, M.M. Gabr, "Discussion of paper by Tong and Lim" J. Roy. Statist. Soc. B , 42 (1980) pp. 278–280
[a17]	T. Teräsvirta, "Specification, estimation, and evaluation of smooth transition autoregressive models" J. Amer. Statist. Assoc. , 89 (1994) pp. 208–218
[a18]	B.Y. Thanoon, "Subset threshold autoregression with applications" J. Time Ser. Anal. , 11 (1990) pp. 75–87
[a19]	D. Tjøsheim, B. Aüestad, "Nonparametric identification of nonlinear time series: selecting significant lags" J. Amer. Statist. Assoc. , 89 (1994) pp. 1410–1419 MR1310231
[a20]	H. Tong, "Some comments on the Canadian lynx data—with discussion" J. Roy. Statist. Soc. A , 140 (1977) pp. 432–435; 448–468
[a21]	H. Tong, "Threshold models in non-linear time series analysis" , Lecture Notes in Statistics , 21 , Springer (1983) Zbl 0527.62083
[a22]	H. Tong, "Non-linear time series: a dynamical system approach" , Clarendon Press (1990) Zbl 0716.62085
[a23]	H. Tong, P. Dabas, "Clusters of time series models: an example" , Techn. Report , Univ. Kent (June, 1989)
[a24]	H. Tong, K.S. Lim, "Threshold autoregression, limit cycles and cyclical data (with discussion)" J. Roy. Statist. Soc. B , 42 (1980) pp. 245–292
[a25]	R.S. Tsay, "Testing and modelling threshold autoregressive processes" J. Amer. Statist. Soc. , 84 (1989) pp. 231–240
[a26]	C. Wong, R. Kohn, "A Bayesian approach to estimating and forecasting additive nonparametric autoregressive models" J. Time Ser. Anal. , 17 (1996) pp. 203–220 MR1381173 Zbl 0845.62068
[a27]	Q. Yao, H. Tong, "Quantifying the influence of initial values on non-linear prediction" J. Roy. Statist. Soc. B , 56 (1994) pp. 701–725 MR1293241
[a28]	Q. Yao, H. Tong, "On subset selection in non-parametric stochastic regression" Statistica Sinica , 4 (1994) pp. 51–70 MR1282865 Zbl 0823.62038
[a29]	G.H. Yu, Y.C. Lin, "A methodology for selecting subset autoregressive time series models" J. Time Ser. Anal. , 12 (1991) pp. 363–373 MR1131008