Threshold autoregressive model
This article is about threshold autoregregressive (TAR) models, which form a specific category in the broader class of regime-switching models. Such models are of interest when you are studying a process that is characterized by recurrent switches between two or more distinct states or regimes that have different means, variances and/or dynamics. For instance, you may believe that the mean and autoregression of a person’s affective well-being change abruptly in response to the level of demands at they experience at their job.
Typical of a TAR model is that the process switches repeatedly between the different regimes, and that this switching process is governed by an observed variable; the latter is referred to as the threshold variable (Tong & Lim, 1980). Each regime itself is characterized by an autoregressive (AR) model. When you have two variables that each serve as the threshold variable for the other variable, this is referred to as a closed-loop threshold autoregressive system (TARSC) model; such models have been used to capture predator-prey cycles (Tong & Lim, 1980), but also for the analysis of dyadic interactions observed in the lab (Hamaker et al., 2009). When the threshold variable is a lagged version of the outcome variable, the process is referred to as a self-exciting threshold autoregressive (SETAR) model; such models have been used in psychological research to captures switches between relapse and recovery in addiction (Warren et al., 2003).
Below you can read more about: 1) the basic TAR model; 2) the SETAR model; 3) several related threshold models; and 4) the challenges involved in estimating and comparing TAR models.
1 The basic TAR model
A TAR model is characterized by a sudden change in one or more aspects of the underlying process whenever the threshold variable passes a specific threshold. In this section you can read more about: a simple TAR model consisting of two regimes with a first-order autoregressive (AR(1)) model for each regime; the long-run means of these regimes; and several other basic TAR models.
1.1 A simple TAR model
In its simplest form, a TAR model has two regimes, and within each regime there is a distinct first-order autoregressive (AR(1)) process. This can be expressed as
\[ y_t = \begin{cases} c_{(1)} + \phi_{(1)} y_{t-1} + \epsilon_t, & \text{if } x_{t-1} < \tau \\[2mm] c_{(2)} + \phi_{(2)} y_{t-1} + \epsilon_t, & \text{if } x_{t-1} \geq \tau \end{cases} \]
where \(x_{t-1}\) is the threshold variable, which determines whether the process is in the first regime or in the second, depending on whether \(x_{t-1}\) is smaller or not than a specific threshold \(\tau\). The difference between the two regimes can be in the intercepts \(c_{(1)}\) versus \(c_{(2)}\), and/or in the autoregression \(\phi_{{(1)}}\) versus \(\phi_{(2)}\). The residual \(\epsilon_t\) comes from a white noise process, which is characterize by a mean of zero and a variance of \(\sigma^2_{\epsilon}\).
Aiko is interested in the regulation of negative affect and how this may be affected by demands at work. She considers using daily diary data in which both negative affect and demands at work are measured.
To model regulation of negative affect, she wants to use a first order autoregressive model, from which she interprets the autoregressive parameter as indicative of the degree to which a person can successfully regulate their negative affect. That is, the closer this parameter is to 0, the less carry-over there is from one occasion to the next, and hence, the better a person can regulate their affect.
Aiko believes that daily fluctuations in job demands do not really have an effect on negative affect, unless it passes a certain threshold; then there may be a sudden increase in average or baseline negative affect and there may be a change in the effectiveness of regulating negative affect. To account for this, Aiko wants to use a TAR model, and she expects that \(c_{(1)} < c_{(2)}\) and \(\phi_{(1)} < \phi_{(2)}\).
1.2 Means of the regimes of a simple TAR model
The long-run means of the regimes may be important features of the model you want to know about: For instance, is it the case that one average negative affect is lower when job demands are low, and is it on average higher when job demands are above a certain threshold. When considering the simple TAR model above, it may be tempting to assume that if \(c_{(1)} < c_{(2)}\), this must imply that the long-run mean in the second regime (when \(x_{t-1} \geq \tau\)) is larger than the mean in the first regime (when \(x_{t-1} < \tau\)). However, this is not necessarily the case.
From the article about autoregressive models, you can see that the long-run mean of a linear AR(1) process is actually a function of the intercept and the autoregressive parameter, that is
\[\mu = \frac{c}{1-\phi}.\]
This shows that the autoregressive parameter also plays a role in getting from the intercept to the long-run mean. But this is only part of the story when considering a TAR model.
The reason for this, is that due to the autoregression, there will be some influence of the features of regime 1 on the long-run features of regime 2, and vice versa. These influences will be relatively minor when: a) the process does not switch that often; b) the autoregressions are not that large; and/or c) the differences in means and variances across the regimes are not that large. But when switches occur frequently and/or the autoregression is substantially, the influence of one regime’s features on that of the other can be quite substantial.
Unfortunately, there are no analytic expression available for the long-run means of a TAR model. If you want to know them, a practical approach is to use the constellation of parameters and simulate a very long TAR process (e.g., of 100 000 occasions), and use this to determine the mean of the observations conditional on the regime. Additionally, you can also determine the variability that characterizes each of the regimes. This can help you to better describe and interpret the regimes from a substantive point of view.
1.3 Other versions of the basic TAR model
The basic TAR model presented above can be extended in various ways. First, there may be more than two regimes, which means there are multiple thresholds. This model can be expressed as
\[ y_t = \begin{cases} c_{(1)} + \phi_{(1)} y_{t-1} + \epsilon_t, & \text{if } x_{t-1} < \tau_1 \\[2mm] c_{(2)} + \phi_{(2)} y_{t-1} + \epsilon_t, & \text{if } \tau_1 \leq x_{t-1} < \tau_2\\[2mm] c_{(3)} + \phi_{(3)} y_{t-1} + \epsilon_t, & \text{if } x_{t-1} \geq \tau_2 \end{cases} \]
This version may be of interest to you when you believe the \(x\) variable can help to distinguish three distinct regimes; for instance, when \(x\) is job demands, you may distinguish between: a regime of very low job demands that may be experienced as boring and demotivating; a regime of healthy job demands, which forms a good match with a person’s abilities and resources; and a regime in which the job demands are very high, which creates frustration, feelings of falling short and failure, and possibly exhaustion. You may also consider more than three regimes, but in practice more regimes will require more data, to be able to estimate such a model.
Second, there can be higher order AR processes in each regime, and the order may differ across the regimes. A general way to express this is
\[ y_t = \begin{cases} c_{(1)} + \phi_{1(1)} y_{t-1} + \phi_{2(1)} y_{t-2} +\dots + \phi_{p_{(1)}(1)} y_{t-p_{(1)}} + \epsilon_t, & \text{if } x_{t-1} < \tau_1 \\[2mm] c_{(2)} + \phi_{1(2)} y_{t-1} + \phi_{2(2)} y_{t-2} +\dots + \phi_{p_{(2)}(2)} y_{t-p_{(2)}} + \epsilon_t, & \text{if } x_{t-1} \geq \tau_2 \end{cases} \]
where \(p_{(1)}\) denotes the order of the AR process in the first regime, and \(p_{(2)}\) denotes the order of the AR process in the second regime. A reason to prefer these higher order AR models in each regime is that this allows for more flexibility in the dynamic patterns over time; you can read more about that in the article about AR models.
Third, the regimes may be characterized by different innovation variances. In this context, it may be important to also realize that if the autoregressive parameters differ across the regimes, this implies that the variance of the outcome variable will also be different, even though the innovation variance remains the same. Yet, to allow for more flexibility, and to avoid bias in the estimation of the autoregressive parameters (Jongerling et al., 2015), you may consider extending the model with regime-specific innovation variances. This can be expresses as
\[ y_t = \begin{cases} c_{(1)} + \phi_{(1)} y_{t-1} + \sigma_{(1)}\epsilon_t, & \text{if } x_{t-1} < \tau \\[2mm] c_{(2)} + \phi_{(2)} y_{t-1} + \sigma_{(2)}\epsilon_t, & \text{if } x_{t-1} \geq \tau \end{cases} \]
where \(\epsilon_t\) has a variance of 1, such that the variance of the residual term \(\sigma_{(1)}\epsilon_t\) is \(\sigma_{(1)}^2\), and the variance of the residual term \(\sigma_{(2)}\epsilon_t\) is \(\sigma_{(2)}^2\). Differences in the innovation variance may be interpreted to mean that the person is for instance more exposed to external factors in one regime than in the other, or that they are more sensitive to such factors in one regime compared to the other.
Aiko believes that when a person is experiencing a period of higher demands at work, this will lead them to become more sensitive to external influences such as minor setbacks or successes throughout the day, support or lack thereof from colleagues, and critical comments and/or compliments from their boss. Instead of measuring all these factors and allowing their effects to differ across the different regimes, Aiko wants to capture this change through including regime-specific innovation variances.
Specifically, Aiko hypothesizes that the second regime, which is entered when there are higher job demands (i.e., \(x_{t-1}>\tau\)), is characterized by a larger innovation variance than the first regime (i.e., \(\sigma_{(1)}^2 < \sigma_{(2)}^2\)).
Fourth, the threshold variable that was considered here was \(x\) at the preceding occasion, that is \(x_{t-1}\). Other lagged versions of \(x\) are also possible, such that the more general expression would be based on having \(x_{t-k}\) as the threshold variable, that is,
\[ y_t = \begin{cases} c_{(1)} + \phi_{(1)} y_{t-1} + \sigma_{(1)}\epsilon_t, & \text{if } x_{t-k} < \tau \\[2mm] c_{(2)} + \phi_{(2)} y_{t-1} + \sigma_{(2)}\epsilon_t, & \text{if } x_{t-k} \geq \tau \end{cases} \]
where \(k\) is referred to as the delay parameter. How much delay there is between a change in the threshold variable and an actual regime switch, may itself be a research question you want to answer.
Finally, you can also consider a lagged version of the outcome as the threshold variable; this version of the TAR model is elaborated on in the next section.
2 Self-exciting threshold autoregressive (SETAR) model
A specific version of the TAR model is when the threshold variable is a lagged version of the outcome variable. The simplest version of this model is
\[ y_t = \begin{cases} c_{(1)} + \phi_{(1)} y_{t-1} + \epsilon_t, & \text{if } y_{t-1} < \tau \\[2mm] c_{(2)} + \phi_{(2)} y_{t-1} + \epsilon_t, & \text{if } y_{t-1} \geq \tau \end{cases} \]
where the switching between the regimes is based on \(y\) at the previous occasion, there are only two regimes, and there is a distinct first-order autoregressive (AR(1)) model for each regime.
When you consider the autoregressive parameter in an AR(1) model as indicative of the degree to which there is a lack of regulation—as an autoregressive parameter closer to 1 implies more carry-over from one occasion to the next, and thus more time is needed to restore equilibrium—then this SETAR model can be used to investigate whether self-regulation is state-dependent. For instance, when you have measured negative affect, you may be interested in investigating whether autoregression is larger—indicating poorer self-regulation—when a certain threshold was passed.
Silvia is interested in the possibility that regulation of affect is state dependent. Specifically, she wants to know whether the autoregression in affect, which has been interpreted as a measure of regulatory weakness, is higher when a person is experiencing negative affect than when they are experiencing positive affect.
To investigate this, Silvia makes use of an existing data set in which the verbal and non-verbal behavior of a person during a conversation with their spouse was represented on a single dimension from extremely negative through neutral to extremely positive affective behavior. Silvia wants to use a self-exciting threshold autoregressive model in which affective behavior at the preceding occasion serves as the threshold variable to distinguish between a positive versus a negative state. The state may differ in terms of their autoregression representing the regulatory weakness, but also the intercept and innovation variance.
With the model Silvia will determine at what value a person switches regimes. Additionally, she expects that the autoregressive parameter in the positive regime is closer to zero than the autoregressive parameter in the negative regime: That would imply that there is more moment-to-moment carry-over when the person is in a negative state than when they are in a positive state.
This example is based on the study by Haan-Rietdijk et al. (2016).
To better understand the various patterns that can be generated with a SETAR model and how the various parameters play a role in this, there are various plots that you can make. These are presented and discussed in more detail below for the very simple SETAR model that was presented above. After that, an interactive tool is presented that allows you to further experiment with different constellations of parameter values of a simple SETAR model.
2.1 Time series plot of a SETAR model
A rather basic plot that you should always make when dealing with time series data, is the plot in which the time series is plotted against time. In Figure 1 you can see an example of this for a dataset that was generated with a SETAR model.
The threshold value is \(\tau=4\), and this is added to the plot as the dashed line. The two regimes have different intercepts, that is, \(c_{(1)}=2\) and \(c_{(2)}=3.5\), while the autoregressive parameter is the same across the two regimes, that is, \(\phi_{(1)} = \phi_{(2)} = 0.4\). Observations from the first regime are shown in purple, and observations from the second regime are shown in green. In addition the sample means per regime are shown as horizontal lines using the same colors. On the right, the densities of the two regimes are plotted.
To see how the parameters affect the patterns, a second example is provided in Figure 2. Here the intercepts are the same across the two regimes, that is \(c_{(1)}=c_{(2)}=2\), but the autoregressive parameter in the second regime is higher than in the first regime, that is \(\phi_{(1)} = 0.2\) while \(\phi_{(2)} = 0.6\). As you can see, the mean in the second regime is again higher than in the first regime, very similar as in the first example. However, this difference now comes from regime differences in the autoregression, not in the intercept.
The examples provided in Figure 1 and Figure 2 quite clearly suggest there are two regimes that the process tends to switch between: Even without the colors, by simply eyeballing the process over time, you may be able to see that it is characterized by two states with different means. However this is not always the case, as illustrated in Figure 3. The process here is similar to the second example in that the intercepts are the same (i.e., \(c_{(1)}=c_{(2)}=2\)), while the autoregressive parameters are different. But this time, the autoregressive parameter in the second regime is lower than in the first regime, that is \(\phi_{(1)} = 0.6\) while \(\phi_{(2)} = 0.2\) (i.e., it is exactly the reverse of the second example).
Here, it is less obvious that there are two distinct regimes in the data generating mechanism. Compared to the previous example, this also shows that how well regimes can be distinguished in terms of their means depends on the combinations of intercepts and autoregressive parameters.
A fourth example is provided in Figure 4, which is a SETAR model in which both the intercept and the autoregression change across regimes. Specifically, the intercept increases from \(c_{(1)}=2\) to \(c_{(2)}=4\), while the autoregressive parameter decreases from \(\phi_{(1)} = 0.6\) to \(\phi_{(2)} = 0.2\) (the latter is the same as in example 3).
Again, it is less obvious from the time series plot that there are two distinct regimes here. This is because these plots are particularly useful to determine whether there are differences in the means of such regimes: If the changes in the parameters are such that the means do not change that much, it can be difficult to see the presence of distinct regimes in these plots. However, there may be other visualizations that are more informative here, as you can see below.
2.2 Autoregressive scatter plot of a SETAR model
A very useful plot when it comes to SETAR models is the autoregressive scatter plot: It is based on plotting \(y_t\) against itself at the preceding occasion, that is, \(y_{t-1}\). In Figure 5 you can see this plot for all four examples that were presented above.
From these scatter plots, you can quite easily see that there is a threshold at \(y_{t-1}=4\): In examples 1-3, the cloud of observations beyond the threshold (i.e., regime 2 in green) is abruptly shifted upward (examples 1 and 2) or downward (example 3) in comparison to the cloud of observation before the threshold (i.e., regime 1 in purple). Only for the fourth example this is less obvious: For this example it may be hard to determine the data come from a SETAR model, because a linear AR(1) model with only one intercept and one autoregressive parameter (i.e., slope) may capture the pattern still pretty well.
Note that, while the time series plot of the third example in Figure 3 provided little reason to assume this process is characterized by distinct regimes, the autoregressive scatter plot provides quite compelling evidence for the existence of two regimes. This is because it focuses explicitly on the relation between \(y_t\) and \(y_{t-1}\), which is a relation that is hard to judge by looking at a time series plot.
Figure 5 also illustrates very clearly what the exact meaning of the regime-specific intercepts and autoregressions are: These are the parameters that determined the regression line when regressing \(y_t\) on \(y_{t-1}\), and indicate at what value this line crosses the y-axis (i.e., the intercept), and what slope characterizes this line (i.e., the autoregression).
2.3 Interactive tool of a SETAR model
You can further explore the kinds of patterns that can be generated with a simple SETAR model with the interactive tool provided below. It is based on a SETAR model with two regimes, each of which is formed by an AR(1) process. The threshold variable is \(y_{t-1}\).
You can specify the intercept, the autoregressive parameter, and the innovation variance for each of the two regimes, and you can change the threshold value. You can also change the number of occasions. Use the refresh button to generate a new series based on the same settings to see how sampling fluctuations may affect the patterns you get to see.
#| '!! shinylive warning !!': |
#| shinylive does not work in self-contained HTML documents.
#| Please set `embed-resources: false` in your metadata.
#| standalone: true
#| viewerHeight: 525
library(shiny)
library(bslib)
## file: app.R
library(shiny)
library(bslib)
library(ggplot2)
library(patchwork)
simulateData <- function(n, c, sigma, phi, tau) {
burnin <- 100
burnin_t <- n + burnin
e1 <- rnorm(burnin_t, sd = sqrt(sigma))
e2 <- rnorm(burnin_t, sd = sqrt(sigma))
y <- numeric(burnin_t)
for (t in 2:burnin_t) {
if (y[t - 1] < tau) {
y[t] <- c[1] + phi[1] * (y[t - 1]) + e1[t]
} else {
y[t] <- c[2] + phi[2] * (y[t - 1]) + e2[t]
}
}
y <- y[(burnin + 1):burnin_t]
dat <- data.frame(y = y, t = 1:n)
return(dat)
}
plotTimeSeries <- function(dat, tau, plot_widths = c(2, 5, 1)) {
c_param <- attr(dat, "c")
phi_param <- attr(dat, "phi")
x_breaks <- pretty(1:nrow(dat))
y_breaks <- pretty(dat$y)
y_lim <- c(min(y_breaks), max(y_breaks))
axis_line_width <- 0.3
dat$y_lag <- c(NA, dat$y[-length(dat$y)])
dat1 <- dat[dat$y_lag < tau, ]
dat2 <- dat[dat$y_lag >= tau, ]
coryy <- round(cor(dat1$y, dat1$y_lag, use = "complete.obs"), 3)
coryy2 <- round(cor(dat2$y, dat2$y_lag, use = "complete.obs"), 3)
meany <- mean(dat$y)
dat$regime <- ifelse(dat$y_lag < tau, "(1)", "(2)")
p_ss <- ggplot(dat, aes(x = y_lag, y = y)) +
geom_vline(xintercept = 0, linewidth = 1, color = "black") +
geom_vline(xintercept = tau, linetype = "dashed", linewidth = 1, color = "grey") +
geom_abline(intercept = c_param[1], slope = phi_param[1], linewidth = 1, colour = "#75216a") +
geom_abline(intercept = c_param[2], slope = phi_param[2], linewidth = 1, colour = "#2fa18f") +
geom_point(aes(color = regime), size = 2.5, shape = 21, fill = "white", stroke = 1.5) +
scale_color_manual(values = c("(1)" = "#75216a", "(2)" = "#2fa18f")) +
scale_y_continuous(breaks = y_breaks) +
scale_x_continuous(breaks = y_breaks) +
coord_cartesian(ylim = y_lim, xlim = y_lim) +
theme_void() +
labs(x = bquote(y[t-1]), y = bquote(y[t])) +
theme(text = element_text(family = "sans", size = 16),
axis.text = element_text(margin = margin(5, 5, 5, 5)),
axis.title = element_text(margin = margin(5, 5, 5, 5)),
axis.ticks = element_line(lineend = "butt", linewidth = axis_line_width),
axis.ticks.length = unit(2.5, "pt"),
panel.spacing = unit(5, units = "pt"),
plot.margin = margin(0, 5, 0, 0),
legend.position = "none") +
geom_segment(x = -Inf, xend = -Inf, y = y_breaks[1],
yend = y_breaks[length(y_breaks)],
linewidth = axis_line_width, lineend = "square") +
geom_segment(y = -Inf, yend = -Inf, x = y_breaks[1],
xend = y_breaks[length(y_breaks)],
linewidth = axis_line_width, lineend = "square") # +
# annotate("text", x = -Inf, y = Inf, hjust = -0.2, vjust = 2,
# label = bquote(hat(rho)[1(1)]~"="~.(coryy)), size = 5) +
# annotate("text", x = -Inf, y = Inf, hjust = -0.2, vjust = 3,
# label = bquote(hat(rho)[1(2)]~"="~.(coryy2)), size = 5)
p_ts <- ggplot(dat) +
#geom_segment(x = x_breaks[1], xend = max(x_breaks),
# y = meany, yend = meany,
# linewidth = 0.3, colour = "gray") +
geom_segment(x = 0, xend = max(dat$t), y = 0, yend = 0,
linewidth = 1, color = "black") +
geom_segment(x = 0, xend = max(dat$t), y = tau, yend = tau,
linetype = "dashed", linewidth = 1, color = "gray") +
geom_line(aes(x = t, y = y)) +
geom_point(aes(x = t, y = y, color = regime), size = 2.5, shape = 21, fill = "white", stroke=1.5, cex=2) +
scale_color_manual(values = c("(1)" = "#75216a", "(2)" = "#2fa18f")) +
scale_y_continuous(breaks = y_breaks) +
scale_x_continuous(breaks = x_breaks) +
coord_cartesian(ylim = y_lim, xlim = c(min(x_breaks), max(x_breaks))) +
theme_void() +
labs(x = "t") +
theme(text = element_text(family = "sans", size = 16),
axis.text.x = element_text(margin = margin(5, 5, 5, 5)),
axis.title.x = element_text(margin = margin(5, 5, 5, 5)),
axis.ticks = element_line(lineend = "butt", linewidth = axis_line_width),
axis.ticks.length = unit(2.5, "pt"),
panel.spacing = unit(5, units = "pt"),
plot.margin = margin(0, 5, 0, 0),
legend.position = "none") +
geom_segment(x = -Inf, xend = -Inf, y = y_breaks[1],
yend = max(y_breaks),
linewidth = axis_line_width, lineend = "square") +
geom_segment(y = -Inf, yend = -Inf, x = x_breaks[1],
xend = max(x_breaks),
linewidth = axis_line_width, lineend = "square")
p_hist <- ggplot(dat) +
geom_histogram(aes(y = y, x = ..density..),
fill = "gray", colour = "black") +
geom_density(aes(y = y)) +
geom_hline(yintercept = meany, colour = "gray", linewidth = 0.5) +
scale_y_continuous(breaks = y_breaks) +
coord_cartesian(ylim = y_lim) +
theme_void() +
labs(x = "", y = "") +
theme(text = element_text(family = "sans", size = 16),
axis.title = element_text(margin = margin(5, 5, 5, 5)),
panel.spacing = unit(5, units = "pt"),
plot.margin = margin(0, 5, 0, 0))
p <- p_ss + p_ts + p_hist + plot_layout(widths = plot_widths)
return(p)
}
ui <- fluidPage(
fluidRow(
column(3, numericInput("n", label = h4("Number of occasions T"), value = 100,
min = 10, max = 1e5, step = 10)),
column(3, numericInput("c1", label = HTML("<h4>Intercept c<sub>(1)</sub></h4>"), value = 1, step = 1)),
column(3, numericInput("var1", label = HTML("<h4>Innovation variance σ<sup>2</sup><sub>ϵ(1)</sub></h4>"),
value = 1, min = 0, max = 100, step = 1)),
column(3, sliderInput("phi1", label = HTML("<h4>Autoregression \U03D5<sub>1(1)</sub></h4>"),
min = -0.999, max = 0.999, value = 0.2, step = 0.05))
),
fluidRow(
column(3, numericInput("tau", label = HTML("<h4>Threshold value τ</h4>"),
value = 3, min = -1e5, max = 1e5, step = 0.1)),
column(3, numericInput("c2", label = HTML("<h4>Intercept c<sub>(2)</sub></h4>"), value = 2, step = 1)),
column(3, numericInput("var2", label = HTML("<h4>Innovation variance σ<sup>2</sup><sub>ϵ(2)</sub></h4>"),
value = 1, min = 0, max = 100, step = 1)),
column(3, sliderInput("phi2", label = HTML("<h4>Autoregression \U03D5<sub>1(2)</sub></h4>"),
min = -0.999, max = 0.999, value = 0.5, step = 0.05))
),
actionButton("refresh", "Refresh"),
plotOutput(outputId = "main_plot", height = "250px")
)
server <- function(input, output) {
dat <- reactive({
input$refresh
d <- simulateData(input$n,
c = c(input$c1, input$c2),
sigma = c(input$var1, input$var2),
phi = c(input$phi1, input$phi2),
tau = input$tau)
attr(d, "c") <- c(input$c1, input$c2)
attr(d, "phi") <- c(input$phi1, input$phi2)
d
})
output$main_plot <- renderPlot({
plotTimeSeries(dat(), tau = input$tau)
})
}
shinyApp(ui = ui, server = server)
The tool shows three visualizations of the data that were generated from the SETAR model with the specified parameter values and number of occasions. The left plot is the autoregressive scatter plot with \(y_t\) plotted against \(y_{t-1}\), using purple for observations from regime 1 and green for observations from regime 2. The regression lines implied by the model parameters are also included using the same coloring. The dashed grey line signifies the threshold.
The middle plot is the time series plotted against time, using the same coloring again. The solid purple line represents the sample mean of observations from the first regime, whereas the solid green line represents the sample mean of observations from the second regime. The dashed grey line represents the threshold value again.
Finally, the right plot shows the distribution as a histogram and a density line, again with a horizontal grey line indicating the observed mean. In contrast to the examples provided above (in Figure 1 to Figure 4), no distinction is made here between observations from the two distinct regimes. Hence, it can be used to see whether the marginal distribution of the empirical data shows evidence for multiple regimes by being bimodal, or that such clues are absent. Bimodality can flag that there are multiple regimes, but as the examples above have also shown, this is not always the case.
2.4 Other versions of the SETAR model
The SETAR model discussed above is rather simple, with only two regimes, each of which is characterized by an AR(1) process. Other versions of the SETAR model include models with more than two regimes, and models with higher order AR processes in the various regimes. Furthermore, the threshold variable may also be a larger lagged version of the outcome, to allow for some delay in the switching process.
Keith, Raymond and Julien have obtained data on daily alcohol consumption from a person who is diagnosed as a heavy drinker. They want to investigate whether a change in alcohol intake is followed by a change in the same direction or in the opposite direction; therefore, they decide to study autoregression in the difference score, rather than in the alcohol consumption itself. Moreover, they believe the dynamics of the process are non-linear, and that these can be captured with a SETAR model with a second-order autoregressive (AR(2)) model in each regime.
Their results indicate that the first regime is characterized by two relatively small negative autoregressive parameters. In contrast, the second regime—which is entered when there was a rather large increase in alcohol consumption—is characterized by an autoregressive parameter close to -1 for the first lag; this implies that there is a strong rebound effect, and any increase is followed by a decrease of almost the same size at the following occasion.
This example is based on Warren et al. (2003).
3 Extensions of the TAR model
There are various ways in which the TAR model has been extended and that may be of interest to you for the study of psychological processes.
First, Tong & Lim (1980) proposed a model in which two variables, \(x\) and \(y\) each are regressed on themselves and each other at the previous occasion(s). Moreover, a lagged version of each variable serves as the threshold variable for the other variable. This model has been called the closed-loop threshold autoregressive system (TARSC), and it has been proposed as a practical approach to study a predator-prey system, which is typically formulated in continuous time using a set of differential equations. The characteristics of such a continuous time system have been shown to be well captured in discrete time with these threshold models (Fan & Yao, 2003).
Hamaker et al. (2009) demonstrated that the nonlinear equations-of-marriage model proposed by Gottman and Murray (Gottman et al., 2002; see also Cook et al., 1995) is actually a TARSC. The equations of marriage were developed to analyze a couple’s interaction in the lab, with a specific interest in identifying the points at which partners abruptly shift their affective behavior in response to each other. For example, one of their models is based on the assumption that a person remains mostly unaffected by the behavior of their spouse when this is fairly neutral; but if their spouse becomes increasingly more negative, at some point this results in a sudden increase in negativity in the person’s behavior as well. Similarly, an increase in positivity does not have an immediate effect, but when the spouse’s positive behavior passes a certain threshold, this lifts the partner’s affective behavior also to a higher level. By showing that these kind of dynamics can be well captured with the TARSC framework, Hamaker et al. (2009) showed that researchers can apply standard TAR-model estimation, model selection and inference methods to study dyadic interaction dynamics.
Second, the TAR model can be extended to allow for autoregressive moving-average (ARMA) models in each regime, resulting in what is known as a TARMA model (Tong, 1993). TARMA models have not been used in psychological research yet, but they should allow for more flexibility in the dynamics (using fewer parameters), just as autoregressive moving-average (ARMA) models allow for more flexibility (with fewer parameters) than pure autoregressive (AR) models.
Third, there are also versions of a TAR model that allow for modeling hysteresis, and that have been referred to as the hysteric TAR model (Li et al., 2015). Hysteresis is the phenomenon that a process switches suddenly when a control variable increases and passes a certain threshold, but that the switch back requires the control variable to decrease further then the value where the initial change occurred. For instance, when job demands increase, there may be a point at which the employee quite suddenly experiences an increase in stress; when the job demands then decrease again, they have to go down much further before the person can experience a sudden decrease in their stress again.
This kind of switching behavior can be modeled within a TAR framework, by allowing the regime the process is in to switch according to the following rules:
\[ S_t = \begin{cases} 1, & \text{if } x_{t-d} < \tau_1 \\[2mm] S_{t-1}, & \text{if } \tau_t \leq x_{t-1} < \tau_2 \\[2mm] 2, & \text{if } x_{t-1} \geq \tau_2 \end{cases} \]
where \(S_t\) indicates whether the process is in the lower or the upper regime (i.e., 1 or 2); \(x_{t-d}\) needs to increase and pass \(\tau_2\) before the system switches from the lower to the upper regime; in contrast, when the system is in the upper regime, \(x_{t-d}\) needs to decrease further than \(\tau_2\), all the way passed \(\tau_1\) before the system switches back to the lower regime. While hysteresis is of interest in various psychological theories, the use of the hysteric TAR model to the study of psychological processes has been limited (de Jong et al., 2023).
4 Estimating and comparing TAR models
Estimating the parameters of a TAR model can be done by reformulating the model so that the estimation problem becomes similar to that of a change point model (although other techniques exist as well, see Fan & Yao, 2003). If the threshold is known, the problem is relatively simple; when the threshold is not known, an exhaustive search is typical to determine where the threshold should be located. This is the same as for a change point model with known or unknown change points.
Comparing models, especially models with different numbers of regimes when the threshold values are unknown and need to be estimated, is a more tricky problem. Again, a similar problem arises in the context of change point models with unknown timing of the change points, and some of the solutions developed there may be used for TAR models as well (Hamaker, 2009).
5 Think more about
TAR models fall within the more extended category of regime-switching models, which also include the change point model, and hidden Markov models (HMMs), such as the regime-switching state-space model (Kim & Nelson, 1999), of which the Markov-switching first-order autoregressive (MS-AR(1)) model may be considered a special case. In comparison to the change point model, the most notable feature of the TAR model is that the threshold variable is not time, but another variable that fluctuates over time. From this perspective, you could think of the change point model as a special version of the TAR model, that is, the version in which the threshold variable is time (i.e., \(x_{t-1} = t\)).
In comparison to the MS-AR model (Hamilton, 1994) and the more general regime-switching state-space model by Kim & Nelson (1999), the most notable feature of the TAR model is that the switching process is not hidden (or latent), but that it is governed by an observed threshold variable. A model that may be a bit in between these two regime-switching models is a version with a HMM, where the switching probabilities depend on an observed variable \(x_t\). Then, the switching is not a deterministic function of \(x_t\) (as it is assumed to be in the TAR model), nor is it completely independent from \(x_t\) (as in the model where the switching probabilities are invariant over time).
It may also be the case that you have obtained direct information about the state a person is in. For instance, you may have data on whether a person is at work or not, or whether they are outside or inside. When you have such a categorical variable about a person’s context, instead of a TAR model, you may consider using an autoregressive moving average model with exogenous predictors (ARMAX model), or dynamic regression analysis (Hyndman & Athanasopoulos, 2021).
While a TAR model is based on switches between distinct regimes and can thereby capture abrupt changes, in general it is a stationary model. This is because, even though the switching occurs over time, it is not triggered by time itself. Moreover, while the switching is reversible, the switches do not occur according to a repetitive deterministic temporal pattern. However, if you use time—or time-of-day or day-of-the-week—as the threshold variable, the model no longer classifies as stationary.
6 Takeaway
A TAR model may be of interest to you when you believe the process you are studying is characterized by two or more distinct states or regimes that the system (e.g., a person or dyad) can be in, and the switches between these regimes can be predicted by an observed variable known as the threshold variable. Particularly interesting versions of a TAR model are the SETAR model in which the outcome variable itself is also the threshold variable, and the TARSC where two variables serve as each others’ threshold variable.
Some typical challenges in TAR modeling concern: a) estimating the threshold values; b) determining how much delay there is between a change in the threshold variable before the system switches regime; c) selecting how many regimes there are; and d) deciding what the order and nature of the processes is in each regime. These challenges may well serve as the actual research questions you want to answer: For instance, your research question may be about how negative a spouse needs to become before a partner responds in anger, or whether a person’s affect regulation can be thought of as switching between two or three distinct states. The problems involved in estimating and selecting TAR models tend to be more challenging than the those encountered when using more basic time series models and require alternative techniques (Fan & Yao, 2003).
When a TAR model has been estimated, you may want to interpret the differences between the regimes in substantive terms; for instance, you may want to be able to describe which of the regimes is characterized by a higher mean, more variability, and/or more inertia (i.e., autoregression). However, it may be challenging to base such statements on the parameters that are obtained, as there tends to be no direct relation between the intercept and the mean, or the variance of the outcome variable and the innovation variance. The analytical expressions that are well known for linear models cannot be applied to the nonlinear TAR model. As a practical solution, you can simulate a long time series using the parameter estimates, and then determine the mean and variance that characterizes each regime empirically.
7 Further reading
We have collected various topics for you to read more about below.
Acknowledgments
This work was supported by the European Research Council (ERC) Consolidator Grant awarded to E. L. Hamaker (ERC-2019-COG-865468).
References
Citation
@article{hamaker2026,
author = {Hamaker, Ellen L. and Berkhout, Sophie W.},
title = {Threshold Autoregressive Model},
journal = {MATILDA},
number = {2026-01-22},
date = {2026-01-22},
url = {https://matilda.fss.uu.nl/articles/tar-model.html},
langid = {en}
}