Learn R Programming

BMTAR

Bayesian Analysis of Multivariate Threshold Autoregressive Models with Missing Data

The R package BMTAR implements parameter estimation using a Bayesian approach for MTAR models (as specific cases, model: AR, VAR and TAR) with missing data using Markov Chain Monte Carlo methods. This package performs the simulation of MTAR process (mtarsim). Estimation of matrix parameters and the threshold values conditional on the autoregressive orders and number of regimes (mtarns). Identification of the autoregressive orders using Bayesian variable selection, together with coefficients and covariance matrices and the threshold values conditional on the number of regimes (mtarstr). Identification of the number of regimes using Metropolised Carlin and Chib or via NAIC criteria (mtarnumreg), to calculate NAIC of any estimated model (mtarNAIC). Estimate missing values together with matrix parameters conditional to threshold values, autoregressive orders and numbers of regimes (mtarmissing). The diagnostic of the residuals in any estimated model can be done (diagnostic_mtar). The package manage several class objects for autoplot and print, functions like (tsregime),(mtaregime) and (mtarinipars) make its construction. Finally, (auto_mtar) its an automatic function that performs all above.

MTAR model

Let and be stochastic processes such that and is a univariate process. follows a MTAR model with threshold variable if:

where is the number of regimes, are the thresholds, which define the regimes. are called output covariates and threshold processes respectively.

Additionally, the innovation process follows a multivariate independent Gaussian zero-mean process with covariance identity matrix it is mutually independent of .

Installation

You can install the development version from Github.

install.packages("devtools")
devtools::install_github("adrincont/BMTAR")

Overview

As mention in the first paragraph lets introduce the objects class and usage in the different functions.

  • tsregime return an object class 'tsregime' which is how the package manage data.
  • mtaregime return an object of class 'regime' use for simulation purposes and as standard presentation of the final estimations.
  • mtarsim return an object of class 'mtarsim' use in autoplot methods. Its practical to conditionate some functions for different known parameters.
  • mtarinipars return an object of class 'regime_inipars' that itself contains an object of class 'tsregime', it is the main object that save known parameters and parameters of the prior distributions for each parameter in a MTAR model. This object needs to be provided in every estimation function.
  • mtarns and mtarstr return an object of class 'regime_model' use in print and autoplot methods, its an standard presentation for estimations done in this functions. It is the object to introduce in mtarNAIC.
  • mtarmissing return an object of class 'regime_missing' for print and autoplot methods.
  • mtarnumreg return an object of class 'regime_number'

Example of use

library(BMTAR)
library(ggplot2)

data(datasim_miss)

data = tsregime(datasim_miss$Yt,datasim_miss$Zt,datasim_miss$Xt)
autoplot.tsregime(data,1)
autoplot.tsregime(data,2)
autoplot.tsregime(data,3)

# Fill in the missing data with the component average
Y_temp = t(datasim_miss$Yt)
meanY = apply(Y_temp,1,mean,na.rm = T)
Y_temp[apply(Y_temp,2,is.na)] = meanY
Y_temp = t(Y_temp)
X_temp = datasim_miss$Xt
meanX = mean(X_temp,na.rm = T)
X_temp[apply(X_temp,2,is.na)] = meanX
Z_temp = datasim_miss$Zt
meanZ = mean(Z_temp,na.rm = T)
Z_temp[apply(Z_temp,2,is.na)] = meanZ

# Estimate the number of regimens with the completed series
data_temp = tsregime(Y_temp,Z_temp,X_temp)
initial = mtarinipars(tsregime_obj = data_temp,list_model = list(l0_max = 3),method = 'KUO')
estim_nr = mtarnumreg(ini_obj = initial,iterprev = 500,niter_m = 500,burn_m = 500, list_m = TRUE,ordersprev = list(maxpj = 2,maxqj = 2,maxdj = 2),parallel = TRUE)
print(estim_nr)

# Estimate the structural and non-structural parameters 
# for the series once we know the number of regimes and some idea of its orders
initial = mtarinipars(tsregime_obj = data_temp,method = 'KUO',
                      list_model = list(pars = list(l = estim_nr$final_m),
                      orders = list(pj = c(2,2))))
estruc = mtarstr(ini_obj = initial,niter = 500,chain = TRUE)
autoplot.regime_model(estruc,1)
autoplot.regime_model(estruc,2)
autoplot.regime_model(estruc,3)
autoplot.regime_model(estruc,4)
autoplot.regime_model(estruc,5)
diagnostic_mtar(estruc)

# With the known structural parameters we estimate the missing data
list_model = list(pars = list(l = estim_nr$final_m,r = estruc$estimates$r[,2],orders = estruc$orders))
initial = mtarinipars(tsregime_obj = datasim_miss,list_model = list_model)
missingest = mtarmissing(ini_obj = initial,chain = TRUE, niter = 500,burn = 500)
print(missingest)
autoplot.regime_missing(missingest,1)
data_c = missingest$tsregim
# ============================================================================================#
# Once the missing data has been estimated, we make the estimates again for all the structural 
# and non-structural parameters.
# ============================================================================================#
initial = mtarinipars(tsregime_obj = data_c,list_model = list(l0_max = 3),method = 'KUO')
estim_nr = mtarnumreg(ini_obj = initial,iterprev = 500,niter_m = 500,burn_m = 500, list_m = TRUE,ordersprev = list(maxpj = 2,maxqj = 2,maxdj = 2))
print(estim_nr)

initial = mtarinipars(tsregime_obj = data_c,method = 'KUO',
list_model = list(pars = list(l = estim_nr$final_m),orders = list(pj = c(2,2))))
estruc = mtarstr(ini_obj = initial,niter = 500,chain = TRUE)
autoplot.regime_model(estruc,1)
autoplot.regime_model(estruc,2)
autoplot.regime_model(estruc,3)
autoplot.regime_model(estruc,4)
autoplot.regime_model(estruc,5)
diagnostic_mtar(estruc)

Other useful examples

MTAR is a general model were it is possible to specificate other kind of models we are familiar with, like

  • Basic auto-regressive model AR(p)
  • Vector auto-regressive model VAR(p)
  • Threshold auto-regressive model TAR(l,pj).
spec/ModelARVARTAR
k1>= 11
Regimes11> 1
Threshold processxx✓

This can be useful when you have missing data in one of this types of models and use BMTAR package for its estimation based on a bayesian approach.

  • AR (with covariates)

If in the MTAR model specification with k = 1, l = 1 and d = 0 we have:

library(mtar)
library(ggplot2)
library(forecast)
# AR = MTAR k = 1, l = 1, Zt = NO
R1 = mtaregime(orders = list(p = 2),Phi = list(phi1 = 0.4,phi2 = 0.3),Sigma = 2)
data = mtarsim(100,list(R1))
ardata = arima.sim(list(ar = c(0.4,0.3),sd = 2),100)
ggpubr::ggarrange(
autoplot(tsregime(ardata)) + ggplot2::labs(title = 'base package'),
autoplot(data$Sim) + ggplot2::labs(title = 'mtar package'),ncol = 2)
arima1 = arima(ts(data$Sim$Yt),c(2,0,0))
parameters = list(l = 1,orders = list(pj = 2))
initial = mtarinipars(tsregime_obj = data$Sim,list_model = list(pars = parameters))
estim1 = mtarns(ini_obj = initial,niter = 1000,chain = TRUE)
print.regime_model(estim1)
ggpubr::ggarrange(
autoplot(estim1,5) + theme(legend.position = 'none') + 
labs(title = 'mtar package'),
ggplot(data = NULL,aes(x = 1:100,y = data$Sim$Yt)) + 
geom_line(col = 'black') + geom_line(data = NULL,
aes(x = 1:100,y = fitted(arima1)),col = "blue") + theme_bw() + 
labs(title = 'forecast package'),ncol = 2)
diagnostic_mtar(estim1)
  • VAR (with covariates)

If in the MTAR model specification with l = 1 and d = 0 we have:

library(mtar)
library(ggplot2)
# VAR = MTAR k > 1, l = 1, Zt = NO
library(vars)
library(BVAR)
library(tsDyn)
R1 = mtaregime(orders = list(p = 1,q = 0,d = 0),
              Phi = list(phi1 = matrix(c(0.3,0.2,0.1,0.4),2,2)),
              Sigma = matrix(c(1,0.5,0.5,1),2,2))
data = mtarsim(100,list(R1))
data2 = tsDyn::VAR.sim(B = matrix(c(0.3,0.2,0.1,0.4),2,2),n = 100,lag = 1,include = c('none'),varcov = matrix(c(1,0.5,0.5,1),2,2))
ggpubr::ggarrange(
autoplot(data$Sim) + labs(title = 'mtar package'),
forecast::autoplot(ts(data2),facets = TRUE) + theme_bw() +
labs(title = 'tsDyn package'),ncol = 2
)
var0 = tsDyn::lineVar(data$Sim$Yt,lag = 1,include = 'none',model = 'VAR')
var1 = vars::VAR(y = data$Sim$Yt,p = 1)
var2 = BVAR::bvar(data = data$Sim$Yt,lags = 1)
parameters = list(l = 1,orders = list(pj = 1))
initial = mtarinipars(tsregime_obj = data$Sim,list_model = list(pars = parameters))
estim1 = mtarns(ini_obj = initial,niter = 1000,chain = TRUE)
estim1$regime
var0
var1$varresult
apply(var2$beta[,,1],2,mean)
apply(var2$beta[,,2],2,mean)
apply(var2$sigma[,,1],2,mean)
apply(var2$sigma[,,2],2,mean)
print.regime_model(estim1)
ggpubr::ggarrange(
autoplot(estim1,5) + theme(legend.position = 'none') + 
labs(title = 'mtar package'),
forecast::autoplot(ts(data$Sim$Yt),facets = TRUE) + theme_bw() +
labs(title = 'tsDyn package') + forecast::autolayer(ts(var0$fitted.values)) +
labs(title = 'tsDyn package') + theme(legend.position = 'none'),ncol = 2)
diagnostic_mtar(estim1)
  • TAR (with covariates)

If in the MTAR model specification with k = 1 we have:

# Example 1, TAR model with 2 regimes
Z = arima.sim(n = 500,list(ar = c(0.5)))
l = 2;r = 0;K = c(2,1)
theta = matrix(c(1,-0.5,0.5,-0.7,-0.3,NA), nrow = l)
H = c(1, 1.5)
X = simu.tar.norm(Z,l,r,K,theta,H)
Yt = tsregim(Yt = X,Zt = Z,r = r)
R1 = mtaregim(orders = list(p = 2),cs = 1,Phi = list(phi1 = -0.5,phi2 = 0.5),
              Sigma = 1)
R2 = mtaregim(orders = list(p = 1),cs = -0.7,Phi = list(phi1 = -0.3),
              Sigma = sqrt(1.5))
YtSim = mtarsim(500,list(R1,R2),r,Zt = Z)
ggpubr::ggarrange(
autoplot(Yt) + ggplot2::labs(title = 'TAR package'),
autoplot(YtSim$Sim) + ggplot2::labs(title = 'mtar package'),ncol = 2)
# number of regimes
res = reg.thr.norm(Z,X)
res$L.est
res$L.prob
res$R.est
res$R.CI
initial = mtarinipars(Yt,list_model = list(l0_min = 2,l0_max = 3),method = 'KUO')
resmtar = mtarnumreg(initial)
# structural parameters
res2 = ARorder.norm(Z,X,l,r)
res2$K.est
res2$K.prob
initial = mtarinipars(Yt,list_model = list(pars = list(l = 2),
orders = list(pj = c(2,2),dj = c(1,1))),method = 'KUO')
res2mtar = mtarstr(initial)
res2mtar$orders
# non-structural parameters
res3 = Param.norm(Z,X,l,r,K) #gibbs
res4 = LS.norm(Z,X,l,r,c(0,0)) #least square
initial = mtarinipars(Yt,list(pars = list(l = 2,orders = list(pj = c(1,1)))))
res3mtar = mtarns(initial)

For more information

You will find the theoretical basis of the method in the documents:

License

This package is free and open source software, licensed under GPL-3.

References

  • Calderon, S. and Nieto, F. (2017) Bayesian analysis of multivariate threshold autoregress models with missing data. Communications in Statistics - Theory and Methods 46 (1):296–318. doi:10.1080/03610926.2014.990758.

Copy Link

Version

Install

install.packages('BMTAR')

Monthly Downloads

336

Version

0.1.1

License

GPL (>= 2)

Maintainer

Andrey Rincon Torres

Last Published

January 19th, 2021

Functions in BMTAR (0.1.1)

autoplot.regime_missing

regime_missing object ggplot for the outputs on the function outputs mtarmissing
datasim_numreg

Multivariate threshold autoregressive process simulation for estimate number of regimes
datasim_miss

Multivariate threshold autoregressive process simulation with missing data
mtarNAIC

Compute NAIC of a MTAR model
lists_ind

Create indicator vector for the regimen of each observation
autoplot.tsregime

tsregime object ggplot for the outputs on the function tsregime
missingest

simulated data
mtaregime

Object class “regime” creation
datasim

Simulated Multivariate threshold autoregressive process simulation
auto_mtar

Estimation of a MTAR model for some data
prodB

Function to make product of elements of a list
print.tsregime

Print tsregime object
dwishartB

Wishart density using Brobdingnag class
repM

Function to create list of matrix objects
hydrodata

Hydrological data of Colombia
mtarns

Estimation of non-structural parameters for MTAR model
mtarnumreg

Estimation of the number of regimes in a MTAR model
tsregime

Creation of class “tsregime” for some data
autoplot

Create a complete ggplot appropriate to a particular data type
print.regime_missing

Print estimates of a regime_missing object of the function output mtarmissing
mtarmissing

Estimation of missing values of observed, covariate and threshold processes
print

print an object appropriate to a particular data type
mtarinipars

Organization and check model specification
print.regime_model

print regime_model object for the function outputs mtarns and mtastr
print.regime_number

print regime_number object for the function outputs mtarnumreg
dmnormB

Multivariate normal density using Brobdingnag class
diagnostic_mtar

Residual diagnosis for model MTAR
mtarsim

Multivariate threshold autoregressive process simulation
mtarstr

Estimation of structural parameters of MTAR model
autoplot.regime_model

regime_model object ggplot for the outputs on the function outputs mtarns and mtastr