Skip to contents

This is a function you can use to simulate time series data for a threshold variable of the HysTAR model. The time series is a (co)sine wave, such that thresholds are crossed in a predictable way. This function is designed to be used in combination with hystar_sim().


  n_t = 100,
  n_switches = 2,
  start_regime = 0,
  start_hyst = FALSE,
  range = c(-1, 1)



The desired length of the simulated time series of z. The actual vector that is returned will contain 10 more time points, see Details. Note that n_t will also be the length of y, when you feed z to hystar_sim.


A scalar indicating the desired number of regime switches. Basically, it is the number of times the variable moves to (and reaches) its minimum or to its maximum. If the thresholds are within the range of z, as they should, this will guarantee the same number of regime switches when the delay parameter of the HysTAR model is greater or equal than the highest order. See Details.


The starting regime of the HysTAR model, 0 (default) or 1.


Logical, should z start in the hysteresis zone? Of course, this also depends on r, and r is not yet specified in this function. Rather, setting start_hyst to TRUE makes z start at in the middle of its range, which makes it easy to set the threshold values "around" the first values of z.


A numeric vector of length 2 indicating the desired range (min, max) of z.


A numeric vector of length n_t. This vector has two attributes "start_regime" and "start_hyst" corresponding to the values you provided. These attributes are used by hystar_sim().


The first value of y that can be predicted in the HysTAR model is at time point \(\max\{d, p\} + 1\), where \(p = \max\{p_0, p_1\}\). This is because we need to observe \(y_{t - p}\) and \(z_{t - d}\). So the first observed value of z that determines a regime is at time point \(\max\{d, p\} + 1 - d\). To make sure that this time point corresponds to the start that you request, z_sim() starts with 10 extra time points. In this way, hystar_sim can select the appropriate time points, based on d and p0, p1.

The HysTAR model

The HysTAR model is defined as:

\( y_t = \begin{cases} \phi_{00} + \phi_{01} y_{t-1} + \cdots + \phi_{0 p_0} y_{t-p_0} + \sigma_{0} \epsilon_{t} \quad \mathrm{if}~R_{t} = 0 \\ \phi_{10} + \phi_{11} y_{t-1} + \cdots + \phi_{1 p_1} y_{t-p_1} + \sigma_{1} \epsilon_{t} \quad \mathrm{if}~R_{t} = 1, \\ \end{cases} \)

with \( R_t = \begin{cases} 0 \quad \quad \mathrm{if} \, z_{t-d} \in (-\infty, r_{0}] \\ R_{t-1} \quad \mathrm{if} \, z_{t-d} \in (r_0, r_1] \\ 1 \quad \quad \mathrm{if} \, z_{t-d} \in (r_1, \infty), \\ \end{cases} \)

where \(p_j\) denotes the order of regime \(j \in \{0,1\}\) with coefficients \(\phi_{j0}, \dots, \phi_{j p_j \in (-1, 1)}\), \(\sigma_{j}\) is the standard deviation of the residuals, and \(d \in \{0, 1, 2, \dots\}\) is a delay parameter. The parameters of primary interest are the thresholds \(r_0 \le r_1\). We let \(t = 0, 1, 2, ..., T\), where \(T\) is the number of observations.


simulated_control_variable <- z_sim()
simulated_hystar_model <- hystar_sim(simulated_control_variable)
fitted_hystar_model <- hystar_fit(simulated_hystar_model$data)