This function simulates multivariate time series generated by a user-specified Threshold Autoregressive (TAR) model.
Usage
simtar(
n,
k = 2,
ars = ars(),
Intercept = TRUE,
trend = c("none", "linear", "quadratic"),
nseason = NULL,
parms,
delay = 0,
thresholds = NULL,
t.series = NULL,
ex.series = NULL,
dist = c("Gaussian", "Student-t", "Hyperbolic", "Laplace", "Slash",
"Contaminated normal", "Skew-Student-t", "Skew-normal"),
skewness = NULL,
extra = NULL,
setar = NULL,
Verbose = TRUE
)Arguments
- n
A positive integer specifying the length of the simulated output series.
- k
A positive integer specifying the dimension of the multivariate output series.
- ars
A list defining the TAR model structure, composed of four elements: the number of regimes (
nregim), the autoregressive order (p), and the maximum lags of the exogenous (q) and threshold (d) series within each regime. This object can be validated using thears()function.- Intercept
An optional logical indicating whether an intercept term is included in each regime.
- trend
An optional character string specifying the degree of deterministic time trend to be included in each regime. Available options are a linear trend (
"linear"), a quadratic trend ("quadratic"), or no trend ("none"). By default,trendis set to"none".- nseason
An optional integer, greater than or equal to 2, specifying the number of seasonal periods. When provided,
nseason - 1seasonal dummy variables are added to the regressors within each regime. By default,nseasonis set toNULL, thereby indicating that the TAR model has no seasonal effects.- parms
A list with one sublist per regime. Each sublist contains two matrices: the first matrix corresponds to the location parameters, and the second matrix corresponds to the scale parameters of the model.
- delay
An optional non-negative integer specifying the delay parameter of the threshold series. By default,
delayis set to0.- thresholds
A numeric vector of length
nregim-1containing the threshold values, sorted in ascending order.- t.series
A matrix containing the values of the threshold series.
- ex.series
A matrix containing the values of the multivariate exogenous series.
- dist
An optional character string specifying the multivariate distribution chosen to model the noise process. Supported options include Gaussian (
"Gaussian"), Student-\(t\) ("Student-t"), Slash ("Slash"), Symmetric Hyperbolic ("Hyperbolic"), Laplace ("Laplace"), and Contaminated Normal ("Contaminated normal"). By default,distis set to"Gaussian".- skewness
An optional numeric vector specifying the skewness parameters of the noise distribution, when applicable.
- extra
An optional value specifying the extra parameter of the noise distribution, when required.
- setar
An optional positive integer indicating which component of the output series is used as the threshold variable. By default,
setaris set toNULL, indicating that the model is not of SETAR type.- Verbose
An optional logical indicating whether a description of the simulated TAR model should be printed. By default,
Verboseis set toTRUE.
Value
A data.frame containing the simulated multivariate output series, and, if
specified, the threshold series and multivariate exogenous series.
References
Vanegas, L.H. and Calderón, S.A. and Rondón, L.M. (2025) Bayesian estimation of a multivariate tar model when the noise process distribution belongs to the class of gaussian variance mixtures. International Journal of Forecasting.
Examples
# \donttest{
###### Simulation of a trivariate TAR model with two regimes
n <- 2000
k <- 3
myars <- ars(nregim=2,p=c(1,2))
Z <- as.matrix(arima.sim(n=n+max(myars$p),list(ar=c(0.5))))
probs <- sort((0.6 + runif(myars$nregim-1)*0.8)*c(1:(myars$nregim-1))/myars$nregim)
dist <- "Student-t"; extra <- 6
parms <- list()
for(j in 1:myars$nregim){
np <- 1 + myars$p[j]*k
parms[[j]] <- list()
parms[[j]]$location <- c(ifelse(runif(np*k)<=0.5,1,-1)*rbeta(np*k,shape1=4,shape2=16))
parms[[j]]$location <- matrix(parms[[j]]$location,np,k)
parms[[j]]$scale <- rgamma(k,shape=1,scale=1)*diag(k)
}
thresholds <- quantile(Z,probs=probs)
out1 <- simtar(n=n, k=k, ars=myars, parms=parms, thresholds=thresholds,
t.series=Z, dist=dist, extra=extra)
#>
#>
#> Sample size :2000 time points
#>
#> Output Series :Y1 | Y2 | Y3
#>
#> Threshold Series :Z with a delay equal to 0
#>
#> Error Distribution :Student-t(6)
#>
#> Number of regimes :2
#>
#> Deterministics :Intercept
#>
#> Autoregressive orders:1, 2
#>
#>
#> Thresholds
#>
#> Regime 1 (-Inf,0.481540646728254]
#> Regime 2 (0.481540646728254,Inf)
#>
#>
#> Regime 1:
#>
#> Autoregressive coefficients
#> Intercept phi_1
#> Y1 -0.1232700 | 0.09574189 0.33472932 0.2218244 |
#> Y2 0.4659198 | 0.31381710 0.07331502 0.4534361 |
#> Y3 -0.0296654 | 0.24993550 -0.11516416 0.1738883 |
#>
#> Scale parameter
#> Y1 Y2 Y3
#> Y1 0.5347256 0.00000000 0.000000000
#> Y2 0.0000000 0.08396186 0.000000000
#> Y3 0.0000000 0.00000000 0.003859399
#>
#> Regime 2:
#>
#> Autoregressive coefficients
#> Intercept phi_1 phi_2
#> Y1 0.2466264 | 0.2394820 0.1451425 -0.1878991 | -0.1731534 -0.2339995
#> Y2 -0.1617798 | -0.1674574 0.2291177 -0.3238317 | 0.1458464 -0.1596195
#> Y3 -0.2219323 | 0.1612170 0.2309067 -0.1729445 | -0.1177893 -0.1121131
#>
#> Y1 -0.2208926 |
#> Y2 0.3060828 |
#> Y3 0.1512029 |
#>
#> Scale parameter
#> Y1 Y2 Y3
#> Y1 0.3407046 0.0000000 0.0000000
#> Y2 0.0000000 0.5927801 0.0000000
#> Y3 0.0000000 0.0000000 0.2615285
str(out1)
#> 'data.frame': 2002 obs. of 5 variables:
#> $ Y1 : num 0.1667 0.0995 1.55 1.5226 0.3912 ...
#> $ Y2 : num 0.671 -0.297 0.929 -0.606 1.184 ...
#> $ Y3 : num -0.248 -0.741 -0.128 0.082 0.41 ...
#> $ Regime: num NA NA 1 2 1 2 1 2 2 2 ...
#> $ Z : num -1.143 -1.142 0.452 0.752 0.369 ...
fit1 <- mtar(~ Y1 + Y2 + Y3 | Z, data=out1, ars=myars, dist=dist,
n.burn=200, n.sim=300, n.thin=2)
summary(fit1)
#>
#>
#> Sample size :1999 time points
#>
#> Output Series :Y1 | Y2 | Y3
#>
#> Threshold Series :Z with a estimated delay equal to 0
#>
#> Error Distribution :Student-t
#>
#> Number of regimes :2
#>
#> Deterministics :Intercept
#>
#> Autoregressive orders:1, 2
#>
#>
#>
#> Thresholds (Mean, HDI.Lower, HDI.Upper)
#>
#> Regime 1 (-Inf,0.48079] (-Inf,0.48065] (-Inf,0.48116]
#> Regime 2 (0.48079,Inf) (0.48065,Inf) (0.48116,Inf)
#>
#>
#> Regime1:
#>
#> Autoregressive coefficients
#> Mean 2(1-PD) HDI.Lower HDI.Upper Mean 2(1-PD)
#> (Intercept) -0.11594 0.00001 -0.17016 -0.05695 | 0.46734 1e-05
#> Y1.lag(1) 0.07352 0.00001 0.02629 0.11821 | 0.31019 1e-05
#> Y2.lag(1) 0.35869 0.00001 0.29871 0.41561 | 0.06723 1e-05
#> Y3.lag(1) 0.20844 0.00667 0.07614 0.31576 | 0.46988 1e-05
#> HDI.Lower HDI.Upper Mean 2(1-PD) HDI.Lower HDI.Upper
#> (Intercept) 0.44981 0.48837 | -0.03054 1e-05 -0.03460 -0.02575
#> Y1.lag(1) 0.29518 0.33086 | 0.25047 1e-05 0.24718 0.25536
#> Y2.lag(1) 0.04053 0.08797 | -0.11405 1e-05 -0.11978 -0.10955
#> Y3.lag(1) 0.43157 0.51862 | 0.16833 1e-05 0.15986 0.17774
#>
#> Scale parameter (Mean, HDI.Lower, HDI.Upper)
#> Y1 Y2 Y3 Y1 Y2 Y3 Y1
#> Y1 0.55078 -0.00527 0.00111 . 0.50039 -0.01770 -0.00127 . 0.60441
#> Y2 -0.00527 0.08344 -0.00072 . -0.01770 0.07643 -0.00177 . 0.00586
#> Y3 0.00111 -0.00072 0.00389 . -0.00127 -0.00177 0.00351 . 0.00389
#> Y2 Y3
#> Y1 0.00586 0.00389
#> Y2 0.09115 0.00026
#> Y3 0.00026 0.00427
#>
#>
#> Regime2:
#>
#> Autoregressive coefficients
#> Mean 2(1-PD) HDI.Lower HDI.Upper Mean 2(1-PD)
#> (Intercept) 0.24814 1e-05 0.18970 0.29835 | -0.07684 0.04667
#> Y1.lag(1) 0.31696 1e-05 0.26373 0.36838 | -0.16504 0.00001
#> Y2.lag(1) 0.11816 1e-05 0.07249 0.16770 | 0.20456 0.00001
#> Y3.lag(1) -0.24397 1e-05 -0.34269 -0.14858 | -0.12809 0.06000
#> Y1.lag(2) -0.18161 1e-05 -0.23281 -0.12060 | 0.19179 0.00001
#> Y2.lag(2) -0.18977 1e-05 -0.26395 -0.13989 | -0.21329 0.00001
#> Y3.lag(2) -0.14105 4e-02 -0.26257 -0.03636 | 0.30165 0.00001
#> HDI.Lower HDI.Upper Mean 2(1-PD) HDI.Lower HDI.Upper
#> (Intercept) -0.14240 0.00342 | -0.23562 1e-05 -0.28703 -0.19243
#> Y1.lag(1) -0.23354 -0.09697 | 0.20154 1e-05 0.15166 0.25037
#> Y2.lag(1) 0.13511 0.26801 | 0.21182 1e-05 0.16190 0.25281
#> Y3.lag(1) -0.24933 0.01367 | -0.14488 1e-05 -0.22876 -0.06692
#> Y1.lag(2) 0.11500 0.26733 | -0.14185 1e-05 -0.18626 -0.09048
#> Y2.lag(2) -0.28229 -0.12325 | -0.12798 1e-05 -0.17715 -0.08046
#> Y3.lag(2) 0.17286 0.43901 | 0.19275 1e-05 0.10528 0.28158
#>
#> Scale parameter (Mean, HDI.Lower, HDI.Upper)
#> Y1 Y2 Y3 Y1 Y2 Y3 Y1
#> Y1 0.32500 0.01540 0.01087 . 0.28719 -0.01325 -0.01272 . 0.37106
#> Y2 0.01540 0.59936 -0.01867 . -0.01325 0.51222 -0.05749 . 0.06359
#> Y3 0.01087 -0.01867 0.25225 . -0.01272 -0.05749 0.22050 . 0.03650
#> Y2 Y3
#> Y1 0.06359 0.03650
#> Y2 0.67959 0.01471
#> Y3 0.01471 0.27969
#>
#>
#> Extra parameter
#> Mean 2(1-PD) HDI.Lower HDI.Upper
#> nu 5.65308 . 4.80652 6.57199
#>
#>
###### Simulation of a trivariate VAR model
n <- 2000
k <- 3
myars <- ars(nregim=1,p=2)
dist <- "Slash"; extra <- 2
parms <- list()
for(j in 1:myars$nregim){
np <- 1 + myars$p[j]*k
parms[[j]] <- list()
parms[[j]]$location <- c(ifelse(runif(np*k)<=0.5,1,-1)*rbeta(np*k,shape1=4,shape2=16))
parms[[j]]$location <- matrix(parms[[j]]$location,np,k)
parms[[j]]$scale <- rgamma(k,shape=1,scale=1)*diag(k)
}
out2 <- simtar(n=n, k=k, ars=myars, parms=parms, dist=dist, extra=extra)
#>
#>
#> Sample size :2000 time points
#>
#> Output Series :Y1 | Y2 | Y3
#>
#> Error Distribution :Slash(2)
#>
#> Number of regimes :1
#>
#> Deterministics :Intercept
#>
#> Autoregressive orders:2 in each regime
#>
#>
#> Regime 1:
#>
#> Autoregressive coefficients
#> Intercept phi_1 phi_2
#> Y1 0.24560209 | 0.2119905 0.07820346 -0.2592281 | -0.2367794 0.2909058
#> Y2 0.30477721 | -0.3630487 0.32088146 0.4332673 | -0.1665745 0.1865442
#> Y3 0.08119541 | 0.1429936 -0.29656954 0.0483543 | 0.1831003 0.1396878
#>
#> Y1 -0.1707754 |
#> Y2 0.3342540 |
#> Y3 0.1439979 |
#>
#> Scale parameter
#> Y1 Y2 Y3
#> Y1 1.533067 0.000000 0.000000
#> Y2 0.000000 1.606902 0.000000
#> Y3 0.000000 0.000000 2.510167
str(out2)
#> 'data.frame': 2002 obs. of 4 variables:
#> $ Y1 : num -0.629 -1.04 1.733 -0.612 -0.852 ...
#> $ Y2 : num 1.848 0.17 -0.926 -1.59 -2.264 ...
#> $ Y3 : num 1.8 1.19 -3.09 -0.56 2.37 ...
#> $ Regime: num NA NA 1 1 1 1 1 1 1 1 ...
fit2 <- mtar(~ Y1 + Y2 + Y3, data=out2, ars=myars, dist=dist,
n.burn=200, n.sim=300, n.thin=2)
summary(fit2)
#>
#>
#> Sample size :2000 time points
#>
#> Output Series :Y1 | Y2 | Y3
#>
#> Error Distribution :Slash
#>
#> Number of regimes :1
#>
#> Deterministics :Intercept
#>
#> Autoregressive orders:2 in each regime
#>
#>
#>
#>
#> Regime1:
#>
#> Autoregressive coefficients
#> Mean 2(1-PD) HDI.Lower HDI.Upper Mean 2(1-PD)
#> (Intercept) 0.31252 1e-05 0.23037 0.40131 | 0.33226 1e-05
#> Y1.lag(1) 0.19444 1e-05 0.17503 0.21354 | -0.35091 1e-05
#> Y2.lag(1) 0.06039 1e-05 0.03862 0.08452 | 0.31873 1e-05
#> Y3.lag(1) -0.25038 1e-05 -0.26396 -0.23198 | 0.43781 1e-05
#> Y1.lag(2) -0.24564 1e-05 -0.26463 -0.22770 | -0.16178 1e-05
#> Y2.lag(2) 0.30538 1e-05 0.28422 0.32488 | 0.19393 1e-05
#> Y3.lag(2) -0.17869 1e-05 -0.19677 -0.15900 | 0.34064 1e-05
#> HDI.Lower HDI.Upper Mean 2(1-PD) HDI.Lower HDI.Upper
#> (Intercept) 0.25235 0.41542 | 0.10327 8e-02 -0.00703 0.21918
#> Y1.lag(1) -0.37301 -0.32950 | 0.15462 1e-05 0.13106 0.18775
#> Y2.lag(1) 0.29541 0.34376 | -0.29456 1e-05 -0.32558 -0.26370
#> Y3.lag(1) 0.42100 0.45713 | 0.07338 1e-05 0.04929 0.09482
#> Y1.lag(2) -0.18155 -0.13918 | 0.18470 1e-05 0.15935 0.21140
#> Y2.lag(2) 0.17173 0.21612 | 0.13520 1e-05 0.10780 0.16113
#> Y3.lag(2) 0.32146 0.35967 | 0.14622 1e-05 0.12153 0.16801
#>
#> Scale parameter (Mean, HDI.Lower, HDI.Upper)
#> Y1 Y2 Y3 Y1 Y2 Y3 Y1 Y2
#> Y1 1.51505 0.00039 0.05663 . 1.37146 -0.07217 -0.01888 . 1.65386 0.07350
#> Y2 0.00039 1.63661 0.03769 . -0.07217 1.46137 -0.06868 . 0.07350 1.77248
#> Y3 0.05663 0.03769 2.64120 . -0.01888 -0.06868 2.38754 . 0.15777 0.13722
#> Y3
#> Y1 0.15777
#> Y2 0.13722
#> Y3 2.84762
#>
#>
#> Extra parameter
#> Mean 2(1-PD) HDI.Lower HDI.Upper
#> nu 1.95811 . 1.83446 2.09005
#>
#>
n <- 5000
k <- 3
myars <- ars(nregim=2,p=c(1,2))
dist <- "Laplace"
parms <- list()
for(j in 1:myars$nregim){
np <- 1 + myars$p[j]*k
parms[[j]] <- list()
parms[[j]]$location <- c(ifelse(runif(np*k)<=0.5,1,-1)*rbeta(np*k,shape1=4,shape2=16))
parms[[j]]$location <- matrix(parms[[j]]$location,np,k)
parms[[j]]$scale <- rgamma(k,shape=1,scale=1)*diag(k)
}
out3 <- simtar(n=n, k=k, ars=myars, parms=parms, delay=2,
thresholds=-1, dist=dist, setar=2)
#>
#>
#> Sample size :5000 time points
#>
#> Output Series :Y1 | Y2 | Y3
#>
#> Threshold Series :Y2 with a delay equal to 2
#>
#> Error Distribution :Laplace
#>
#> Number of regimes :2
#>
#> Deterministics :Intercept
#>
#> Autoregressive orders:1, 2
#>
#>
#> Thresholds
#>
#> Regime 1 (-Inf,-1]
#> Regime 2 (-1,Inf)
#>
#>
#> Regime 1:
#>
#> Autoregressive coefficients
#> Intercept phi_1
#> Y1 0.2092999 | 0.1152458 0.1765656 0.2234638 |
#> Y2 -0.3013037 | -0.1671852 0.1751160 0.2663608 |
#> Y3 0.3356962 | -0.1486716 -0.1166652 0.1083313 |
#>
#> Scale parameter
#> Y1 Y2 Y3
#> Y1 0.7763125 0.000000 0.0000000
#> Y2 0.0000000 1.112193 0.0000000
#> Y3 0.0000000 0.000000 0.4155577
#>
#> Regime 2:
#>
#> Autoregressive coefficients
#> Intercept phi_1 phi_2
#> Y1 -0.25393002 | 0.08585486 -0.20140673 -0.1701490 | 0.2396091 -0.1431742
#> Y2 -0.30620616 | -0.08977061 -0.08581429 -0.1533585 | -0.1711748 0.2461271
#> Y3 -0.09011262 | 0.22673250 0.13620635 -0.2570248 | -0.2668896 -0.2861827
#>
#> Y1 0.19863693 |
#> Y2 0.09365131 |
#> Y3 0.12983483 |
#>
#> Scale parameter
#> Y1 Y2 Y3
#> Y1 1.121163 0.000000 0.000000
#> Y2 0.000000 1.204285 0.000000
#> Y3 0.000000 0.000000 0.812355
str(out3)
#> 'data.frame': 5002 obs. of 4 variables:
#> $ Y1 : num 0.337 -0.855 2.906 1.948 1.32 ...
#> $ Y2 : num 0.7762 0.7069 -2.3025 -1.0387 0.0735 ...
#> $ Y3 : num 0.917 1.494 0.632 -7.253 -6.44 ...
#> $ Regime: num NA NA 2 2 1 1 2 1 1 1 ...
fit3 <- mtar(~ Y1 + Y2 + Y3, data=out3, ars=myars, dist=dist,
n.burn=200, n.sim=300, n.thin=2, setar=2)
summary(fit3)
#>
#>
#> Sample size :4999 time points
#>
#> Output Series :Y1 | Y2 | Y3
#>
#> Threshold Series :Y2 with a estimated delay equal to 2
#>
#> Error Distribution :Laplace
#>
#> Number of regimes :2
#>
#> Deterministics :Intercept
#>
#> Autoregressive orders:1, 2
#>
#>
#>
#> Thresholds (Mean, HDI.Lower, HDI.Upper)
#>
#> Regime 1 (-Inf,-1.00199] (-Inf,-1.00578] (-Inf,-0.99975]
#> Regime 2 (-1.00199,Inf) (-1.00578,Inf) (-0.99975,Inf)
#>
#>
#> Regime1:
#>
#> Autoregressive coefficients
#> Mean 2(1-PD) HDI.Lower HDI.Upper Mean 2(1-PD)
#> (Intercept) 0.19502 1e-05 0.13631 0.25080 | -0.25243 1e-05
#> Y1.lag(1) 0.11072 1e-05 0.08978 0.13411 | -0.18131 1e-05
#> Y2.lag(1) 0.17555 1e-05 0.15385 0.19775 | 0.18126 1e-05
#> Y3.lag(1) 0.22753 1e-05 0.20262 0.25530 | 0.24860 1e-05
#> HDI.Lower HDI.Upper Mean 2(1-PD) HDI.Lower HDI.Upper
#> (Intercept) -0.32600 -0.18510 | 0.35770 1e-05 0.31744 0.39419
#> Y1.lag(1) -0.20910 -0.15658 | -0.15021 1e-05 -0.16684 -0.13064
#> Y2.lag(1) 0.15951 0.20273 | -0.10732 1e-05 -0.12054 -0.09356
#> Y3.lag(1) 0.21283 0.27729 | 0.10241 1e-05 0.07993 0.11933
#>
#> Scale parameter (Mean, HDI.Lower, HDI.Upper)
#> Y1 Y2 Y3 Y1 Y2 Y3 Y1
#> Y1 0.75276 0.00405 -0.01387 . 0.68325 -0.04480 -0.04388 . 0.81022
#> Y2 0.00405 1.09746 0.01039 . -0.04480 1.00336 -0.01935 . 0.05453
#> Y3 -0.01387 0.01039 0.43260 . -0.04388 -0.01935 0.39644 . 0.01612
#> Y2 Y3
#> Y1 0.05453 0.01612
#> Y2 1.19977 0.04881
#> Y3 0.04881 0.47394
#>
#>
#> Regime2:
#>
#> Autoregressive coefficients
#> Mean 2(1-PD) HDI.Lower HDI.Upper Mean 2(1-PD)
#> (Intercept) -0.31017 1e-05 -0.36417 -0.25457 | -0.33121 1e-05
#> Y1.lag(1) 0.08169 1e-05 0.06491 0.09560 | -0.09527 1e-05
#> Y2.lag(1) -0.20473 1e-05 -0.21885 -0.19076 | -0.08064 1e-05
#> Y3.lag(1) -0.17980 1e-05 -0.19812 -0.16397 | -0.14697 1e-05
#> Y1.lag(2) 0.23979 1e-05 0.22452 0.25634 | -0.16228 1e-05
#> Y2.lag(2) -0.12380 1e-05 -0.14772 -0.10082 | 0.23546 1e-05
#> Y3.lag(2) 0.19915 1e-05 0.18173 0.21854 | 0.08742 1e-05
#> HDI.Lower HDI.Upper Mean 2(1-PD) HDI.Lower HDI.Upper
#> (Intercept) -0.38323 -0.27017 | -0.09090 1e-05 -0.13928 -0.03801
#> Y1.lag(1) -0.11293 -0.07915 | 0.22506 1e-05 0.21349 0.23736
#> Y2.lag(1) -0.09500 -0.06432 | 0.13478 1e-05 0.12341 0.14716
#> Y3.lag(1) -0.16467 -0.12813 | -0.26546 1e-05 -0.28144 -0.24950
#> Y1.lag(2) -0.18314 -0.14399 | -0.26527 1e-05 -0.27897 -0.25242
#> Y2.lag(2) 0.20599 0.26056 | -0.28183 1e-05 -0.30170 -0.26570
#> Y3.lag(2) 0.06748 0.10712 | 0.13200 1e-05 0.11783 0.14594
#>
#> Scale parameter (Mean, HDI.Lower, HDI.Upper)
#> Y1 Y2 Y3 Y1 Y2 Y3 Y1
#> Y1 1.08795 -0.00528 0.00740 . 1.01306 -0.04202 -0.03004 . 1.15693
#> Y2 -0.00528 1.18071 -0.01342 . -0.04202 1.10958 -0.05433 . 0.03906
#> Y3 0.00740 -0.01342 0.78171 . -0.03004 -0.05433 0.72741 . 0.04014
#> Y2 Y3
#> Y1 0.03906 0.04014
#> Y2 1.26718 0.02363
#> Y3 0.02363 0.83617
#>
#>
# }