Consider the lynx data available in the \(\mathrm{R}\) package tsDyn. These data consist of annual totals of Canadian lynx trapped in the Mackenzie River district of NW Canada during 1821-1934. Define \(y_{t}\) as a \(\log 10\) transformation of the original series. The regression model to be considered is

\[y_{t}=\beta_{0}+\beta_{1} y_{t-1}+\beta_{2} y_{t-1}+u_{t}, \quad u_{t} \sim N\left(0, \sigma^{2}\right)\]

where the regression coefficients and the residual variance are subject to different regimes, and SETAR threshold functions are defined as

\[K_{t}=y_{t-2}-\delta_{\mathrm{j}}, \quad j=1, \ldots, J-1\]

Compare models with regimes involving \(J=2\) and \(J=3\) states, for example using the LPML, and obtain the proportion of periods spent in each regime. A code for the two regime case (conditioning on \(y_{1}\) and \(y_{2}\), and without a stationarity constraint on the AR parameters) is

model for (t in 3:114) (y[t] ~ dnorm (mu [t], tau [S[t]]) LL[t]