Introduction to the mtarm Package"

In mtarm: Bayesian Estimation of Multivariate Threshold Autoregressive Models

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 6,
  fig.height = 4
)

Multivariate Threshold Autoregressive (TAR) models

The mtarm package provides a computational tool designed for Bayesian estimation, inference, and forecasting in multivariate Threshold Autoregressive (TAR) models. These models provide a versatile approach for modeling nonlinear multivariate time series and include multivariate Self-Exciting Threshold Autoregressive (SETAR) and Vector Autoregressive (VAR) models as particular cases [@VANEGAS2025]. The package accommodates a broad class of innovation distributions beyond the Gaussian assumption, such as Student-$t$, slash, symmetric hyperbolic, Laplace, contaminated normal, skew-normal, and skew-$t$ distributions, thereby enabling robust modeling of heavy tails, asymmetry, and other non-Gaussian characteristics.

Installation

Install from GitHub

remotes::install_github("lhvanegasp/mtarm")

Install from CRAN

install.packages("mtarm")

Application: Temperature, precipitation, and two river flows in Iceland

Dataset

The data are available in the object `iceland.rf` and were obtained from [@T90], who provided a detailed description of the geographical and meteorological characteristics of the rivers and analyzed each series individually. Subsequently, [@T98] conducted a bivariate analysis of the same dataset. The focus is on the bivariate time series ${(Y_{1,t},Y_{2,t})^{\top}}{t\geq 1}$, where $Y{1,t}$ and $Y_{2,t}$ denote the daily river flow (in cubic meters per second, ${m}^3/{s}$) of the Jökulsá Eystri and Vatnsdalsá rivers, respectively. The sample covers the period from 1972 to 1974, comprising 1095 observations. The exogenous variables include daily precipitation $X_t$, measured in millimeters (${mm}$), and temperature $Z_t$, measured in degrees Celsius ($^\circ\mathrm{C}$), both recorded at the meteorological station in Hveravellir. Precipitation corresponds to the accumulated rainfall from 9:00 A.M. of the previous day to 9:00 A.M. of the current day.

library(mtarm)
data(iceland.rf)       
str(iceland.rf)     

summary(iceland.rf[,-5])

plot(ts(as.matrix(iceland.rf[,-5])), main="Iceland")

Model specification

Following [@T98], the series are modeled using a $\mathrm{TAR}(2; p=(15,15), q=(4,4), d=(2,2))$ specification given by

$$ Y_t= \sum_{j=1}^{2} I!\left(Z_{t-h}\in(c_{j-1},c_j]\right) \left(! \phi_0^{^{(j)}} +\sum_{i=1}^{15}\boldsymbol{\phi}i^{^{(j)}}Y{t-i} +\sum_{i=1}^{4}\boldsymbol{\beta}i^{^{(j)}}X{t-i} +\sum_{i=1}^{2}\delta_i^{^{(j)}}Z_{t-i} +\epsilon_t^{^{(j)}} !\right) $$

where $\epsilon_t^{^{(j)}}$ is the error term. The last 55 observations (from November 7 to December 31, 1974), corresponding to $5\%$ of the sample, are excluded from the estimation stage and reserved for out-of-sample forecast evaluation. The following code requests the estimation for the $\mathrm{TAR}(2; p=(15,15), q=(4,4), d=(2,2))$ specification under Gaussian, Student-$t$, and Laplace error distributions.

Parameter estimation

set.seed(09102)
fits <- mtar_grid(~ Jokulsa + Vatnsdalsa | Temperature | Precipitation,
                  data=iceland.rf, subset={Date<="1974-11-06"},                           
                  row.names=Date, nregim.min=2, nregim.max=2, p.min=15,                 
                  p.max=15, q.min=4, q.max=4, d.min=2, d.max=2,                           
                  n.burnin=500, n.sim=400, n.thin=2, ssvs=TRUE,
                  dist=c("Gaussian","Student-t","Laplace"),
                  plan_strategy="multisession")

fits

Model comparison using forecast accuracy measures

Adjusted within-sample

The following code requests Deviance Information Criterion (DIC) [@spiegelhalter02; @spiegelhalter14] and Watanabe-Akaike Information Criterion (WAIC) [@w10] values.

DIC(fits)

WAIC(fits)

Out-of-sample with standard $h$-step-ahead forecasting

In addition, the following code provides the median of the log-score [@Good1952], the Energy Score (ES) [@GneitingEtAl2008TEST]—a multivariate extension of the Continuous Ranked Probability Score (CRPS)[@MathesonWinkler1976; @GrimitGneitingBerrocalJohnson2006]—and the Absolute Percentage Error (APE), all computed from the observed and forecasted values for the last 55 observations.

newdata <- subset(iceland.rf, Date>"1974-11-06") 
set.seed(09102)

oos <- out_of_sample(fits, newdata=newdata, n.ahead=nrow(newdata), FUN=median) 
oos[,c(1,2,5,6)]

Out-of-sample with rolling-origin forecasting and fixed parameters

set.seed(09102)
oos2 <- out_of_sample(fits, newdata=newdata, n.ahead=nrow(newdata), 
                      rolling=5, FUN=median) 

for(i in 1:length(oos2)){                   
   cat("\n",i,"-step-ahead\n")   
   print(oos2[[i]][,c(1,2,5,6)]) 
}                            

Summary of the best model

summary(fits[["Laplace.2.15.4.2"]])

Residuals

res <- residuals(fits[["Laplace.2.15.4.2"]])  

par(mfrow=c(1,2)) 
qqnorm(res[["full"]], pch=20, col="blue", main="") 
abline(0, 1, lty=3) 
hist(res[["full"]], freq=FALSE, xlab="Quantile-type residual",
     ylab="Density", main="") 
curve(dnorm(x), col="blue", add=TRUE)

par(mfrow=c(1,2))  
acf(res[["full"]], col="blue", main="")
pacf(res[["full"]], col="blue", main="")

Forecasting

pred <- predict(fits[["Laplace.2.15.4.2"]], newdata=newdata,
                n.ahead=nrow(newdata), row.names=Date, credible=0.8)

head(pred$summary)
tail(pred$summary)

Summary statistics

fitmcmc <- coda::as.mcmc(fits[["Laplace.2.15.4.2"]])
summary(fitmcmc)

Convergence diagnostics

Geweke statistic

geweke_diagTAR(fits[["Laplace.2.15.4.2"]])

Effective sample size

effectiveSize_TAR(fits[["Laplace.2.15.4.2"]]) 

mtarm documentation built on June 12, 2026, 5:07 p.m.