Title: | Analysis of Nonlinear Time Series Through Threshold Autoregressive Moving Average Models (TARMA) Models |
---|---|
Description: | Routines for nonlinear time series analysis based on Threshold Autoregressive Moving Average (TARMA) models. It provides functions and methods for: TARMA model fitting and forecasting, including robust estimators, see Goracci et al. JBES (2025) <doi:10.1080/07350015.2024.2412011>; tests for threshold effects, see Giannerini et al. JoE (2024) <doi:10.1016/j.jeconom.2023.01.004>, Goracci et al. Statistica Sinica (2023) <doi:10.5705/ss.202021.0120>, Angelini et al. (2024) <doi:10.48550/arXiv.2308.00444>; unit-root tests based on TARMA models, see Chan et al. Statistica Sinica (2024) <doi:10.5705/ss.202022.0125>. |
Authors: | Simone Giannerini [aut, cre] , Greta Goracci [aut] |
Maintainer: | Simone Giannerini <[email protected]> |
License: | GPL (>= 2) |
Version: | 0.5-1 |
Built: | 2024-12-08 07:17:52 UTC |
Source: | CRAN |
The data is taken from Table 1 of (Andrews 2003), which provides asymptotic critical values for sup Wald, LM, and LR tests for parameter instability. Critical values are given for degrees of freedom \(p=1,\dots,20\). They can be used as asymptotic critical values for the following tests:
Provided pb = 1- pa
.
ACValues
ACValues
ACValues
A matrix with 13 rows and 62 columns. The first two columns contain the parameters, whereas the remaining 60 columns contain 3 critical values (at level 10%, 5%, and 1%) for each value of the parameter \(p\):
pi
\(\pi_0\) in the paper, corresponds to pa
.
lambda
\(\lambda\) in the paper, not relevant here.
1-90
\(p=1\), critical level 10%.
1-95
\(p=1\), critical level 5%.
1-99
\(p=1\), critical level 1%.
...
20-90
\(p=20\), critical level 10%.
20-95
\(p=20\), critical level 5%.
20-99
\(p=20\), critical level 1%.
Note that \(p\) are the degrees of freedom and correspond to the total number of tested parameter in the above tests.
Andrews DWK (2003). “Tests for Parameter Instability and Structural Change with Unknown Change Point: A Corrigendum.” Econometrica, 71(1), 395-397. doi:10.1111/1468-0262.00405.
Plots a time series model fit, possibly together with its forecast and forecast bands.
## S3 method for class 'tsfit' plot( x, fit, plot.fit = TRUE, fore = NULL, lcols = list(series = "lightblue", fit = 4, pred = "red4", band = "red"), ptype = list(series = 20, fit = 1, pred = 16, band = 20), ltype = list(series = 1, fit = 1, pred = 1, band = 1), lwdt = 2, ... )
## S3 method for class 'tsfit' plot( x, fit, plot.fit = TRUE, fore = NULL, lcols = list(series = "lightblue", fit = 4, pred = "red4", band = "red"), ptype = list(series = 20, fit = 1, pred = 16, band = 20), ltype = list(series = 1, fit = 1, pred = 1, band = 1), lwdt = 2, ... )
x |
A time series |
fit |
A time series model fitted upon |
plot.fit |
Logical. If |
fore |
Forecast derived from |
lcols |
List of colours for the plots. |
ptype |
List of point types ( |
ltype |
List of line types ( |
lwdt |
A common line width for the plots. |
... |
Additional graphical parameters. |
No return value, called for side effects
Simone Giannerini, [email protected]
Greta Goracci, [email protected]
TARMA.fit
and TARMA.fit2
for TARMA modelling.
predict.TARMA
for prediction and forecasting.
## a TARMA(1,1,1,1) model set.seed(13) x <- TARMA.sim(n=200, phi1=c(0.5,-0.5), phi2=c(0.0,0.5), theta1=-0.5, theta2=0.7, d=1, thd=0.2) fit1 <- TARMA.fit(x,tar1.lags = 1, tar2.lags = 1, tma1.lags = 1, tma2.lags = 1, d=1, threshold=0.2) xp1 <- predict(fit1,x,n.ahead=5) # plots both the fitted and the forecast plot.tsfit(x,fit=fit1,fore=xp1);grid(); # plots only the forecast plot.tsfit(x,fit=fit1,plot.fit=FALSE,fore=xp1);grid(); # plots only the fitted plot.tsfit(x,fit=fit1);grid();
## a TARMA(1,1,1,1) model set.seed(13) x <- TARMA.sim(n=200, phi1=c(0.5,-0.5), phi2=c(0.0,0.5), theta1=-0.5, theta2=0.7, d=1, thd=0.2) fit1 <- TARMA.fit(x,tar1.lags = 1, tar2.lags = 1, tma1.lags = 1, tma2.lags = 1, d=1, threshold=0.2) xp1 <- predict(fit1,x,n.ahead=5) # plots both the fitted and the forecast plot.tsfit(x,fit=fit1,fore=xp1);grid(); # plots only the forecast plot.tsfit(x,fit=fit1,plot.fit=FALSE,fore=xp1);grid(); # plots only the fitted plot.tsfit(x,fit=fit1);grid();
Forecasting with TARMA models
## S3 method for class 'TARMA' predict( object, x, n.ahead = 0, n.sim = 1000, quant = c(0.05, 0.95), pred.matrix = FALSE, ... )
## S3 method for class 'TARMA' predict( object, x, n.ahead = 0, n.sim = 1000, quant = c(0.05, 0.95), pred.matrix = FALSE, ... )
object |
A |
x |
The fitted time series. |
n.ahead |
The number of steps ahead for which prediction is required. |
n.sim |
The number of Monte Carlo replications used to simulate the prediction density. |
quant |
Vector of quantiles (in the interval |
pred.matrix |
Logical. if |
... |
Additional arguments. |
If n.ahead = 0
it gives the fitted values from the model.
If the fit is from TARMA.fit2
and includes covariates, these are ignored.
A list with components pred.matrix
, pred
, and pred.interval
. The latter two are ts
objects that contain the prediction and the quantiles of the prediction density, respectively.
If pred.matrix = TRUE
then the prediction density from which the quantiles are computed is also returned.
Simone Giannerini, [email protected]
Greta Goracci, [email protected]
Giannerini S, Goracci G (2021). “Estimating and Forecasting with TARMA models.” University of Bologna.
TARMA.fit
and TARMA.fit2
for TARMA modelling. plot.tsfit
for plotting TARMA fits and forecasts.
## a TARMA(1,1,1,1) model set.seed(13) x1 <- TARMA.sim(n=200, phi1=c(0.5,-0.5), phi2=c(0.0,0.5), theta1=-0.5, theta2=0.7, d=1, thd=0.2) fit1 <- TARMA.fit(x1, method='L-BFGS-B',tar1.lags = 1, tar2.lags = 1, tma1.lags = 1, tma2.lags = 1, d=1, threshold=0.2) xp1 <- predict(fit1,x1,n.ahead=2) xp1
## a TARMA(1,1,1,1) model set.seed(13) x1 <- TARMA.sim(n=200, phi1=c(0.5,-0.5), phi2=c(0.0,0.5), theta1=-0.5, theta2=0.7, d=1, thd=0.2) fit1 <- TARMA.fit(x1, method='L-BFGS-B',tar1.lags = 1, tar2.lags = 1, tma1.lags = 1, tma2.lags = 1, d=1, threshold=0.2) xp1 <- predict(fit1,x1,n.ahead=2) xp1
Methods for TARMA fits
## S3 method for class 'TARMA' print(x, digits = max(3L, getOption("digits") - 3L), se = TRUE, ...) ## S3 method for class 'TARMA' coef(object, ...) ## S3 method for class 'TARMA' vcov(object, ...) ## S3 method for class 'TARMA' residuals(object, ...)
## S3 method for class 'TARMA' print(x, digits = max(3L, getOption("digits") - 3L), se = TRUE, ...) ## S3 method for class 'TARMA' coef(object, ...) ## S3 method for class 'TARMA' vcov(object, ...) ## S3 method for class 'TARMA' residuals(object, ...)
x |
A |
digits |
Number of decimal digits for the output. |
se |
Logical. if |
... |
Further parameters. |
object |
A |
No return value, called for side effects
TARMA.fit
and TARMA.fit2
for TARMA modelling, plot.tsfit
for plotting TARMA fits and forecasts.
Methods for TARMA tests
## S3 method for class 'TARMAtest' print(x, ...)
## S3 method for class 'TARMAtest' print(x, ...)
x |
A |
... |
Further parameters passed to |
Print method for TARMAtest
objects. Prints the results using the method for class htest
and adds critical values extracted from ACValues
for the test for threshold nonlinearity
and supLMQur
for the unit root test against the TARMA alternative.
Note that the bootstrap version of the tests also print the bootstrap p-value.
No return value, called for side effects
Simone Giannerini, [email protected]
Greta Goracci, [email protected]
The data provides asymptotic null critical values for the unit root supLM test described in (Chan et al. 2024). They are used with the following tests:
Provided pb = 1- pa
.
supLMQur
supLMQur
supLMQur
A 4-dimensional array that contains 4 critical values (at level 10%, 5%, 1%, 0.1%) for each combination of
pa
Lower bound for the threshold range. From 0.01 to 0.4
th
MA(1) parameter.
n
Sample size.
Chan K, Giannerini S, Goracci G, Tong H (2024). “Testing for threshold regulation in presence of measurement error.” Statistica Sinica, 34(3), 1413-1434. doi:10.5705/ss.202022.0125, https://doi.org/10.5705/ss.202022.0125.
Heteroskedasticity robust supremum Lagrange Multiplier test for a AR specification versus a TAR specification. Includes the classic (non robust) AR versus TAR test.
TAR.test(x, pa = 0.25, pb = 0.75, ar.ord, d = 1)
TAR.test(x, pa = 0.25, pb = 0.75, ar.ord, d = 1)
x |
A univariate time series. |
pa |
Real number in |
pb |
Real number in |
ar.ord |
Order of the AR part. |
d |
Delay parameter. Defaults to |
Implements a heteroskedasticity robust asymptotic supremum Lagrange Multiplier test to test an AR specification versus a TAR specification.
This is an asymptotic test and the value of the test statistic has to be compared with the critical values tabulated
in (Goracci et al. 2021) or (Andrews 2003).
Both the non-robust supLM
and the robust supLMh
statistics are returned.
An object of class TARMAtest
with components:
statistic
A named vector with the values of the classic supLM
and robust supLMh
statistics.
parameter
A named vector: threshold
is the value that maximises the Lagrange Multiplier values.
test.v
Matrix of values of the LM statistic for each threshold given in thd.range
.
thd.range
Range of values of the threshold.
fit
The null model: AR fit over x
.
sigma2
Estimated innovation variance from the AR fit.
data.name
A character string giving the name of the data.
prop
Proportion of values of the series that fall in the lower regime.
p.value
The p-value of the test. It is NULL
for the asymptotic test.
method
A character string indicating the type of test performed.
d
The delay parameter.
pa
Lower threshold quantile.
dfree
Effective degrees of freedom. It is the number of tested parameters.
Simone Giannerini, [email protected]
Greta Goracci, [email protected]
Goracci G, Giannerini S, Chan K, Tong H (2023). “Testing for threshold effects in the TARMA framework.” Statistica Sinica, 33(3), 1879-1901. doi:10.5705/ss.202021.0120.
Andrews DWK (2003). “Tests for Parameter Instability and Structural Change with Unknown Change Point: A Corrigendum.” Econometrica, 71(1), 395-397. doi:10.1111/1468-0262.00405.
TAR.test.B
for the bootstrap version of the test.
TARMA.test
for the (robust) ARMA vs TARMA asymptotic version of the test, which includes also the AR vs TAR test, with different defaults.
TARMAGARCH.test
for the robust version of the ARMA vs TARMA test that assumes GARCH innovations.
TARMA.sim
to simulate from a TARMA process.
set.seed(123) ## a TAR(1,1) --------------- x1 <- TARMA.sim(n=100, phi1=c(0.5,-0.5), phi2=c(0.0,0.8), theta1=0, theta2=0, d=1, thd=0.2) TAR.test(x1, ar.ord=1, d=1) ## a AR(1) ---------------- x2 <- arima.sim(n=100, model=list(order=c(1,0,0), ar=0.5)) TAR.test(x2, ar.ord=1, d=1)
set.seed(123) ## a TAR(1,1) --------------- x1 <- TARMA.sim(n=100, phi1=c(0.5,-0.5), phi2=c(0.0,0.8), theta1=0, theta2=0, d=1, thd=0.2) TAR.test(x1, ar.ord=1, d=1) ## a AR(1) ---------------- x2 <- arima.sim(n=100, model=list(order=c(1,0,0), ar=0.5)) TAR.test(x2, ar.ord=1, d=1)
Implements various bootstrap supremum Lagrange Multiplier tests for a AR specification versus a TAR specification.
TAR.test.B( x, B = 1000, pa = 0.25, pb = 0.75, ar.ord, d = 1, btype = c("iid", "wb.h", "wb.r", "wb.n"), ... )
TAR.test.B( x, B = 1000, pa = 0.25, pb = 0.75, ar.ord, d = 1, btype = c("iid", "wb.h", "wb.r", "wb.n"), ... )
x |
A univariate time series. |
B |
Integer. Number of bootstrap resamples. Defaults to 1000. |
pa |
Real number in |
pb |
Real number in |
ar.ord |
Order of the AR part. |
d |
Delay parameter. Defaults to |
btype |
Bootstrap type, can be one of |
... |
Additional arguments to be passed to |
Implements the bootstrap version of TAR.test
the supremum Lagrange Multiplier test to test an AR specification versus a TARMA specification.
The option btype
specifies the type of bootstrap as follows:
iid
Residual iid bootstrap. See (Giannerini et al. 2022), (Giannerini et al. 2024).
wb.h
Stochastic permutation of (Hansen 1996).
wb.r
Residual wild bootstrap with Rademacher auxiliary distribution. See (Giannerini et al. 2022), (Giannerini et al. 2024).
wb.n
Residual wild bootstrap with Normal auxiliary distribution. See (Giannerini et al. 2022), (Giannerini et al. 2024).
A list of class htest
with components:
statistic
The value of the supLM statistic.
parameter
A named vector: threshold
is the value that maximises the Lagrange Multiplier values.
test.v
Vector of values of the LM statistic for each threshold given in thd.range
.
thd.range
Range of values of the threshold.
fit
The null model: AR fit over x
.
sigma2
Estimated innovation variance from the AR fit.
data.name
A character string giving the name of the data.
prop
Proportion of values of the series that fall in the lower regime.
p.value
The bootstrap p-value of the test.
method
A character string indicating the type of test performed.
Tb
The bootstrap null distribution.
Simone Giannerini, [email protected]
Greta Goracci, [email protected]
Giannerini S, Goracci G, Rahbek A (2022). “The validity of bootstrap testing in the threshold framework.” doi:10.48550/ARXIV.2201.00028, https://arxiv.org/abs/2201.00028.
Giannerini S, Goracci G, Rahbek A (2024). “The validity of bootstrap testing in the threshold framework.” Journal of Econometrics, 239(1), 105379. ISSN 0304-4076, doi:10.1016/j.jeconom.2023.01.004, Climate Econometrics.
Goracci G, Giannerini S, Chan K, Tong H (2023). “Testing for threshold effects in the TARMA framework.” Statistica Sinica, 33(3), 1879-1901. doi:10.5705/ss.202021.0120.
Giannerini S, Goracci G (2021). “Estimating and Forecasting with TARMA models.” University of Bologna.
Hansen BE (1996). “Inference When a Nuisance Parameter Is Not Identified Under the Null Hypothesis.” Econometrica, 64(2), 413–430. ISSN 00129682, 14680262, https://doi.org/10.2307/2171789.
TAR.test
for the heteroskedastic robust asymptotic test. TARMAGARCH.test
for the
robust version of the test with respect to GARCH innovations. TARMA.sim
to simulate from a TARMA process.
## a TAR(1,1) where the threshold effect is on the AR parameters set.seed(123) x1 <- TARMA.sim(n=100, phi1=c(0.5,-0.5), phi2=c(0.0,0.8), theta1=0, theta2=0, d=1, thd=0.2) TAR.test.B(x1, ar.ord=1, d=1) TAR.test.B(x1, ar.ord=1, d=1, btype='wb.r') TAR.test.B(x1, ar.ord=1, d=1, btype='wb.h') ## a AR(1) x2 <- arima.sim(n=100, model=list(order = c(1,0,0),ar=0.5)) TAR.test.B(x2, ar.ord=1, d=1) TAR.test.B(x2, ar.ord=1, d=1, btype='wb.r') TAR.test.B(x2, ar.ord=1, d=1, btype='wb.h')
## a TAR(1,1) where the threshold effect is on the AR parameters set.seed(123) x1 <- TARMA.sim(n=100, phi1=c(0.5,-0.5), phi2=c(0.0,0.8), theta1=0, theta2=0, d=1, thd=0.2) TAR.test.B(x1, ar.ord=1, d=1) TAR.test.B(x1, ar.ord=1, d=1, btype='wb.r') TAR.test.B(x1, ar.ord=1, d=1, btype='wb.h') ## a AR(1) x2 <- arima.sim(n=100, model=list(order = c(1,0,0),ar=0.5)) TAR.test.B(x2, ar.ord=1, d=1) TAR.test.B(x2, ar.ord=1, d=1, btype='wb.r') TAR.test.B(x2, ar.ord=1, d=1, btype='wb.h')
Implements a Least Squares fit of full subset two-regime TARMA(p1,p2,q1,q2)
model to a univariate time series
TARMA.fit( x, tar1.lags = c(1), tar2.lags = c(1), tma1.lags = c(1), tma2.lags = c(1), threshold = NULL, d = 1, pa = 0.25, pb = 0.75, method = c("L-BFGS-B", "solnp", "lbfgsb3c", "robust", "trimmed"), alpha = 0, qu = c(0.05, 0.95), innov = c("norm", "student"), optim.control = list(trace = 0), irls.control = list(maxiter = 100, tol = 1e-04), ... )
TARMA.fit( x, tar1.lags = c(1), tar2.lags = c(1), tma1.lags = c(1), tma2.lags = c(1), threshold = NULL, d = 1, pa = 0.25, pb = 0.75, method = c("L-BFGS-B", "solnp", "lbfgsb3c", "robust", "trimmed"), alpha = 0, qu = c(0.05, 0.95), innov = c("norm", "student"), optim.control = list(trace = 0), irls.control = list(maxiter = 100, tol = 1e-04), ... )
x |
A univariate time series. |
tar1.lags |
Vector of AR lags for the lower regime. It can be a subset of |
tar2.lags |
Vector of AR lags for the upper regime. It can be a subset of |
tma1.lags |
Vector of MA lags for the lower regime. It can be a subset of |
tma2.lags |
Vector of MA lags for the upper regime. It can be a subset of |
threshold |
Threshold parameter. If |
d |
Delay parameter. Defaults to |
pa |
Real number in |
pb |
Real number in |
method |
Optimization/fitting method, can be one of |
alpha |
Real positive number. Tuning parameter for robust estimation. Only used if |
qu |
Quantiles for (initial) trimmed estimation. Tuning parameter for robust estimation. Only used if |
innov |
Innovation density for robust estimation. can be one of |
optim.control |
List of control parameters for the main optimization method. |
irls.control |
List of control parameters for the irls optimization method (see details). |
... |
Additional arguments. |
Implements the Least Squares fit of the following two-regime TARMA(p1,p2,q1,q2)
process:
\[X_{t} = \left\lbrace
\begin{array}{ll}
\phi_{1,0} + \sum_{i \in I_1} \phi_{1,i} X_{t-i} + \sum_{j \in M_1} \theta_{1,j} \varepsilon_{t-j} + \varepsilon_{t} & \mathrm{if } X_{t-d} \leq \mathrm{thd} \\
&\\
\phi_{2,0} + \sum_{i \in I_2} \phi_{2,i} X_{t-i} + \sum_{j \in M_2} \theta_{2,j} \varepsilon_{t-j} + \varepsilon_{t} & \mathrm{if } X_{t-d} > \mathrm{thd}
\end{array}
\right. \]
where \(\phi_{1,i}\) and \(\phi_{2,i}\) are the TAR parameters for the lower and upper regime, respectively, and
I1 = tar1.lags
and I2 = tar2.lags
are the corresponding vectors of TAR lags.
\(\theta_{1,j}\) and \(\theta_{2,j}\) are the TMA parameters
and \(j \in M_1, M_2\), where M1 = tma1.lags
and M2 = tma2.lags
, are the vectors of TMA lags.
The most demanding routines have been reimplemented in Fortran and dynamically loaded.
A list of class TARMA
with components:
fit
- List with the following components
coef
- Vector of estimated parameters which can be extracted by the coef method.
sigma2
- Estimated innovation variance.
var.coef
- The estimated variance matrix of the coefficients coef, which can be extracted by the vcov method
residuals
- Vector of residuals from the fit.
nobs
- Effective sample size used for fitting the model.
se
- Standard errors for the parameters. Note that they are computed conditionally upon the threshold so that they are generally smaller than the true ones.
thd
- Estimated threshold.
aic
- Value of the AIC for the minimised least squares criterion over the threshold range.
bic
- Value of the BIC for the minimised least squares criterion over the threshold range.
rss
- Minimised value of the target function. Coincides with the residual sum of squares for ordinary least squares estimation.
rss.v
- Vector of values of the rss over the threshold range.
thd.range
- Vector of values of the threshold range.
d
- Delay parameter.
phi1
- Estimated AR parameters for the lower regime.
phi2
- Estimated AR parameters for the upper regime.
theta1
- Estimated MA parameters for the lower regime.
theta2
- Estimated MA parameters for the upper regime.
tlag1
- TAR lags for the lower regime
tlag2
- TAR lags for the upper regime
mlag1
- TMA lags for the lower regime
mlag2
- TMA lags for the upper regime
method
- Estimation method.
innov
- Innovation density model.
alpha
- Tuning parameter for robust estimation.
qu
- Tuning parameter for robust estimation.
call
- The matched call.
convergence
- Convergence code from the optimization routine.
innovpar
- Parameter vector for the innovation density. Defaults to NULL
.
method
has the following options:
L-BFGS-B
Calls the corresponding method of optim
. Linear ergodicity constraints are imposed.
solnp
Calls the function solnp
. It is a nonlinear optimization using augmented Lagrange method with
linear and nonlinear inequality bounds. This allows to impose all the ergodicity constraints so that in theory it always
return an ergodic solution. In practice the solution should be checked since this is a local solver and there is no guarantee
that the minimum has been reached.
lbfgsb3c
Calls the function lbfgsb3c
in package lbfgsb3c
. Improved version of the L-BFGS-B in optim
.
robust
Robust M-estimator of Ferrari and La Vecchia (Ferrari and La-Vecchia 2011). Based on the L-BFGS-B in optim
and an additional iterative re-weighted least squares step to estimate the robust weights.
Uses the tuning parameters alpha
and qu
. Robust standard errors are derived from the sandwich estimator of the variance/covariance matrix of the estimates.
The IRLS step can be controlled through the parameters maxiter
(maximum number of iterations) and tol
(target tolerance). These can be passed using irls.control
.
trimmed
Experimental: Estimator based on trimming the sample using the tuning parameters qu
(lower and upper quantile).
Where possible, the conditions for ergodicity and invertibility are imposed to the optimization routines but there is no guarantee that the solution will be ergodic and invertible so that it is advisable to check the fitted parameters.
Simone Giannerini, [email protected]
Greta Goracci, [email protected]
Giannerini S, Goracci G (2021). “Estimating and Forecasting with TARMA models.” University of Bologna.
Chan K, Goracci G (2019). “On the Ergodicity of First-Order Threshold Autoregressive Moving-Average Processes.” J. Time Series Anal., 40(2), 256-264.
Goracci G, Ferrari D, Giannerini S, Ravazzolo F (2023). “Robust estimation for Threshold Autoregressive Moving-Average models.” Free University of Bolzano, University of Bologna. doi:10.48550/ARXIV.2211.08205.
Ferrari D, La-Vecchia D (2011). “On robust estimation via pseudo-additive information.” Biometrika, 99(1), 238-244. ISSN 0006-3444, doi:10.1093/biomet/asr061.
TARMA.fit2
for Maximum Likelihood estimation of TARMA
models with common MA part. print.TARMA
for print methods for TARMA
fits.
predict.TARMA
for prediction and forecasting. plot.tsfit
for plotting TARMA fits and forecasts.
## a TARMA(1,1,1,1) model set.seed(13) x <- TARMA.sim(n=200, phi1=c(0.5,-0.5), phi2=c(0.0,0.5), theta1=-0.5, theta2=0.7, d=1, thd=0.2) fit1 <- TARMA.fit(x,tar1.lags=1, tar2.lags=1, tma1.lags=1, tma2.lags=1, d=1) ## -------------------------------------------------------------------------- ## In the following examples the threshold is fixed to speed up computations ## -------------------------------------------------------------------------- ## -------------------------------------------------------------------------- ## Least Squares fit ## -------------------------------------------------------------------------- set.seed(26) n <- 200 y <- TARMA.sim(n=n, phi1=c(0.6,0.6), phi2=c(-1.0,0.4), theta1=-0.7, theta2=0.5, d=1, thd=0.2) fit1 <- TARMA.fit(y,tar1.lags=1, tar2.lags=1, tma1.lags=1, tma2.lags=1, d=1, threshold=0.2) fit1 ## --------------------------------------------------------------------------- ## Contaminate the data with one additive outlier ## --------------------------------------------------------------------------- x <- y # contaminated series x[54] <- x[54] + 10 ## --------------------------------------------------------------------------- ## Compare the non-robust LS fit with the robust fit ## --------------------------------------------------------------------------- fitls <- TARMA.fit(x,tar1.lags=1, tar2.lags=1, tma1.lags=1, tma2.lags=1, d=1, threshold=0.2) fitrob <- TARMA.fit(x,tar1.lags=1, tar2.lags=1, tma1.lags=1, tma2.lags=1, d=1, method='robust',alpha=0.7,qu=c(0.1,0.95), threshold=0.2) par.true <- c(0.6,0.6,-1,0.4,-0.7,0.5) pnames <- c("int.1", "ar1.1", "int.2", "ar2.1", "ma1.1", "ma2.1") names(par.true) <- pnames par.ls <- round(fitls$fit$coef,2) # Least Squares par.rob <- round(fitrob$fit$coef,2) # robust rbind(par.true,par.ls,par.rob)
## a TARMA(1,1,1,1) model set.seed(13) x <- TARMA.sim(n=200, phi1=c(0.5,-0.5), phi2=c(0.0,0.5), theta1=-0.5, theta2=0.7, d=1, thd=0.2) fit1 <- TARMA.fit(x,tar1.lags=1, tar2.lags=1, tma1.lags=1, tma2.lags=1, d=1) ## -------------------------------------------------------------------------- ## In the following examples the threshold is fixed to speed up computations ## -------------------------------------------------------------------------- ## -------------------------------------------------------------------------- ## Least Squares fit ## -------------------------------------------------------------------------- set.seed(26) n <- 200 y <- TARMA.sim(n=n, phi1=c(0.6,0.6), phi2=c(-1.0,0.4), theta1=-0.7, theta2=0.5, d=1, thd=0.2) fit1 <- TARMA.fit(y,tar1.lags=1, tar2.lags=1, tma1.lags=1, tma2.lags=1, d=1, threshold=0.2) fit1 ## --------------------------------------------------------------------------- ## Contaminate the data with one additive outlier ## --------------------------------------------------------------------------- x <- y # contaminated series x[54] <- x[54] + 10 ## --------------------------------------------------------------------------- ## Compare the non-robust LS fit with the robust fit ## --------------------------------------------------------------------------- fitls <- TARMA.fit(x,tar1.lags=1, tar2.lags=1, tma1.lags=1, tma2.lags=1, d=1, threshold=0.2) fitrob <- TARMA.fit(x,tar1.lags=1, tar2.lags=1, tma1.lags=1, tma2.lags=1, d=1, method='robust',alpha=0.7,qu=c(0.1,0.95), threshold=0.2) par.true <- c(0.6,0.6,-1,0.4,-0.7,0.5) pnames <- c("int.1", "ar1.1", "int.2", "ar2.1", "ma1.1", "ma2.1") names(par.true) <- pnames par.ls <- round(fitls$fit$coef,2) # Least Squares par.rob <- round(fitrob$fit$coef,2) # robust rbind(par.true,par.ls,par.rob)
Maximum Likelihood fit of a two-regime TARMA(p1,p2,q,q)
model with common MA parameters, possible common AR parameters and possible covariates.
TARMA.fit2( x, ar.lags = NULL, tar1.lags = c(1), tar2.lags = c(1), ma.ord = 1, sma.ord = 0L, period = NA, threshold = NULL, d = 1, pa = 0.25, pb = 0.75, thd.var = NULL, include.int = TRUE, x.reg = NULL, optim.control = list(), ... )
TARMA.fit2( x, ar.lags = NULL, tar1.lags = c(1), tar2.lags = c(1), ma.ord = 1, sma.ord = 0L, period = NA, threshold = NULL, d = 1, pa = 0.25, pb = 0.75, thd.var = NULL, include.int = TRUE, x.reg = NULL, optim.control = list(), ... )
x |
A univariate time series. |
ar.lags |
Vector of common AR lags. Defaults to |
tar1.lags |
Vector of AR lags for the lower regime. It can be a subset of |
tar2.lags |
Vector of AR lags for the upper regime. It can be a subset of |
ma.ord |
Order of the MA part (also called |
sma.ord |
Order of the seasonal MA part (also called |
period |
Period of the seasonal MA part (also called |
threshold |
Threshold parameter. If |
d |
Delay parameter. Defaults to |
pa |
Real number in |
pb |
Real number in |
thd.var |
Optional exogenous threshold variable. If |
include.int |
Logical. If |
x.reg |
Covariates to be included in the model. These are passed to |
optim.control |
List of control parameters for the optimization method. |
... |
Additional arguments passed to |
Fits the following two-regime TARMA
process with optional components: linear AR
part, seasonal MA
and covariates.
\[X_{t} = \phi_{0} + \sum_{h \in I} \phi_{h} X_{t-h} + \sum_{l=1}^Q \Theta_{l} \epsilon_{t-ls} + \sum_{j=1}^q \theta_{j} \epsilon_{t-j} + \sum_{k=1}^K \delta_{k} Z_{k} + \epsilon_{t} + \left\lbrace
\begin{array}{ll}
\phi_{1,0} + \sum_{i \in I_1} \phi_{1,i} X_{t-i} & \mathrm{if } X_{t-d} \leq \mathrm{thd} \\
&\\
\phi_{2,0} + \sum_{i \in I_2} \phi_{2,i} X_{t-i} & \mathrm{if } X_{t-d} > \mathrm{thd}
\end{array}
\right. \]
where \(\phi_h\) are the common AR parameters and \(h\) ranges in I = ar.lags
. \(\theta_j\) are the common MA parameters and \(j = 1,\dots,q\)
(q = ma.ord
), \(\Theta_l\) are the common seasonal MA parameters and \(l = 1,\dots,Q\) (Q = sma.ord
)
\(\delta_k\) are the parameters for the covariates. Finally, \(\phi_{1,i}\) and \(\phi_{2,i}\) are the TAR parameters
for the lower and upper regime, respectively and I1 = tar1.lags
I2 = tar2.lags
are the vector of TAR lags.
A list of class TARMA
with components:
fit
- The output of the fit. It is a arima
object.
aic
- Value of the AIC for the minimised least squares criterion over the threshold range.
bic
- Value of the BIC for the minimised least squares criterion over the threshold range.
aic.v
- Vector of values of the AIC over the threshold range.
thd.range
- Vector of values of the threshold range.
d
- Delay parameter.
thd
- Estimated threshold.
phi1
- Estimated AR parameters for the lower regime.
phi2
- Estimated AR parameters for the upper regime.
theta1
- Estimated MA parameters for the lower regime.
theta2
- Estimated MA parameters for the upper regime.
delta
- Estimated parameters for the covariates x.reg
.
tlag1
- TAR lags for the lower regime
tlag2
- TAR lags for the upper regime
mlag1
- TMA lags for the lower regime
mlag2
- TMA lags for the upper regime
arlag
- Same as the input slot ar.lags
include.int
- Same as the input slot include.int
se
- Standard errors for the parameters. Note that they are computed conditionally upon the threshold so that they are generally smaller than the true ones.
rss
- Minimised residual sum of squares.
method
- Estimation method.
call
- The matched call.
Simone Giannerini, [email protected]
Greta Goracci, [email protected]
Giannerini S, Goracci G (2021). “Estimating and Forecasting with TARMA models.” University of Bologna.
Chan K, Goracci G (2019). “On the Ergodicity of First-Order Threshold Autoregressive Moving-Average Processes.” J. Time Series Anal., 40(2), 256-264.
TARMA.fit
for Least Square estimation of full subset TARMA
models. print.TARMA
for print methods for TARMA
fits.
predict.TARMA
for prediction and forecasting.
## a TARMA(1,1,1,1) set.seed(127) x <- TARMA.sim(n=100, phi1=c(0.5,-0.5), phi2=c(0,0.8), theta1=0.5, theta2=0.5, d=1, thd=0.2) fit1 <- TARMA.fit2(x, tar1.lags=1, tar2.lags=1, ma.ord=1, d=1) ## Showcase of the fit with covariates --- ## simulates from a TARMA(3,3,1,1) model with common MA parameter ## and common AR(1) and AR(2) parameters. Only the lag 3 parameter varies across regimes set.seed(212) n <- 300 x <- TARMA.sim(n=n, phi1=c(0.5,0.3,0.2,0.4), phi2=c(0.5,0.3,0.2,-0.2), theta1=0.4, theta2=0.4, d=1, thd=0.2, s1=1, s2=1) ## FIT 1: estimates lags 1,2,3 as threshold lags --- fit1 <- TARMA.fit2(x, ma.ord=1, tar1.lags=c(1,2,3), tar2.lags=c(1,2,3), d=1) ## FIT 2: estimates lags 1 and 2 as fixed AR and lag 3 as the threshold lag fit2 <- TARMA.fit2(x, ma.ord=1, tar1.lags=c(3), tar2.lags=c(3), ar.lags=c(1,2), d=1) ## FIT 3: creates lag 1 and 2 and fits them as covariates --- z1 <- lag(x,-1) z2 <- lag(x,-2) fit3 <- TARMA.fit2(x, ma.ord=1, tar1.lags=c(3), tar2.lags=c(3), x.reg=ts.intersect(z1,z2), d=1) ## FIT 4: estimates lag 1 as a covariate, lag 2 as fixed AR and lag 3 as the threshold lag fit4 <- TARMA.fit2(x, ma.ord = 1, tar1.lags=c(3), tar2.lags=c(3), x.reg=z1, ar.lags=2, d=1)
## a TARMA(1,1,1,1) set.seed(127) x <- TARMA.sim(n=100, phi1=c(0.5,-0.5), phi2=c(0,0.8), theta1=0.5, theta2=0.5, d=1, thd=0.2) fit1 <- TARMA.fit2(x, tar1.lags=1, tar2.lags=1, ma.ord=1, d=1) ## Showcase of the fit with covariates --- ## simulates from a TARMA(3,3,1,1) model with common MA parameter ## and common AR(1) and AR(2) parameters. Only the lag 3 parameter varies across regimes set.seed(212) n <- 300 x <- TARMA.sim(n=n, phi1=c(0.5,0.3,0.2,0.4), phi2=c(0.5,0.3,0.2,-0.2), theta1=0.4, theta2=0.4, d=1, thd=0.2, s1=1, s2=1) ## FIT 1: estimates lags 1,2,3 as threshold lags --- fit1 <- TARMA.fit2(x, ma.ord=1, tar1.lags=c(1,2,3), tar2.lags=c(1,2,3), d=1) ## FIT 2: estimates lags 1 and 2 as fixed AR and lag 3 as the threshold lag fit2 <- TARMA.fit2(x, ma.ord=1, tar1.lags=c(3), tar2.lags=c(3), ar.lags=c(1,2), d=1) ## FIT 3: creates lag 1 and 2 and fits them as covariates --- z1 <- lag(x,-1) z2 <- lag(x,-2) fit3 <- TARMA.fit2(x, ma.ord=1, tar1.lags=c(3), tar2.lags=c(3), x.reg=ts.intersect(z1,z2), d=1) ## FIT 4: estimates lag 1 as a covariate, lag 2 as fixed AR and lag 3 as the threshold lag fit4 <- TARMA.fit2(x, ma.ord = 1, tar1.lags=c(3), tar2.lags=c(3), x.reg=z1, ar.lags=2, d=1)
TARMA(p1,p2,q1,q2)
processSimulates from the following two-regime TARMA(p1,p2,q1,q2)
process:
TARMA.sim( n, phi1, phi2, theta1, theta2, d = 1, thd = 0, s1 = 1, s2 = 1, rand.gen = rnorm, innov = rand.gen(n, ...), n.start = 500, xstart, start.innov = rand.gen(n.start, ...), ... )
TARMA.sim( n, phi1, phi2, theta1, theta2, d = 1, thd = 0, s1 = 1, s2 = 1, rand.gen = rnorm, innov = rand.gen(n, ...), n.start = 500, xstart, start.innov = rand.gen(n.start, ...), ... )
n |
Length of the series. |
phi1 |
Vector of |
phi2 |
Vector of |
theta1 |
Vector of |
theta2 |
Vector of |
d |
Delay parameter. Defaults to |
thd |
Threshold parameter. Defaults to |
s1 |
Innovation variance for the lower regime. Defaults to |
s2 |
Innovation variance for the upper regime. Defaults to |
rand.gen |
Optional: a function to generate the innovations. Defaults to |
innov |
Optional: a time series of innovations. If not provided, |
n.start |
Length of the burn-in period. Defaults to |
xstart |
Initial condition as a named list: |
start.innov |
Optional: a time series of innovations for the burn-in period. |
... |
Additional arguments for |
Note that the parameters are not checked for ergodicity.
A time series object of class ts
generated from the above model.
Simone Giannerini, [email protected]
Greta Goracci, [email protected]
Giannerini S, Goracci G (2021). “Estimating and Forecasting with TARMA models.” University of Bologna.
## a TARMA(1,1,1,1) model set.seed(123) x <- TARMA.sim(n=100, phi1=c(0.5,-0.5), phi2=c(0.0,0.8), theta1=-0.5, theta2=0.5, d=1, thd=0.2) ## a TARMA(1,2,1,1) model x <- TARMA.sim(n=100,phi1=c(0.5,-0.5,0),phi2=c(0,0.5,0.3),theta1=-0.5,theta2=0.5,d=1,thd=0.2)
## a TARMA(1,1,1,1) model set.seed(123) x <- TARMA.sim(n=100, phi1=c(0.5,-0.5), phi2=c(0.0,0.8), theta1=-0.5, theta2=0.5, d=1, thd=0.2) ## a TARMA(1,2,1,1) model x <- TARMA.sim(n=100,phi1=c(0.5,-0.5,0),phi2=c(0,0.5,0.3),theta1=-0.5,theta2=0.5,d=1,thd=0.2)
Heteroskedasticity robust supremum Lagrange Multiplier tests for a ARMA specification versus a TARMA specification. Includes the AR versus TAR test.
TARMA.test( x, pa = 0.25, pb = 0.75, ar.ord, ma.ord, ma.fixed = TRUE, d, thd.range, method = "CSS-ML", ... )
TARMA.test( x, pa = 0.25, pb = 0.75, ar.ord, ma.ord, ma.fixed = TRUE, d, thd.range, method = "CSS-ML", ... )
x |
A univariate time series, either a |
pa |
Real number in |
pb |
Real number in |
ar.ord |
Order of the AR part. |
ma.ord |
Order of the MA part. |
ma.fixed |
Logical. Only applies to testing ARMA vs TARMA. If |
d |
Delay parameter. Defaults to |
thd.range |
Vector of optional user defined threshold range. If missing then |
method |
Fitting method to be passed to |
... |
Additional arguments to be passed to |
Implements an asymptotic supremum Lagrange Multiplier test to test an ARMA specification versus a TARMA specification.
Both the non-robust supLM
and the robust supLMh
statistics are returned.
If ma.fixed=TRUE
(the default), the AR parameters are tested whereas the MA parameters are fixed. If ma.fixed=FALSE
both the AR and the MA parameters are tested.
This is an asymptotic test and the value of the test statistic has to be compared with the critical values tabulated in (Goracci et al. 2021) and (Andrews 2003). These are automatically computed and printed by print.TARMAtest
.
If ma.ord=0
then the AR versus TAR test is used. Note that when method='CSS'
, this is equivalent to TAR.test
, which uses least squares.
An object of class TARMAtest
with components:
statistic
The value of the supLM
statistic and its robust version supLMh
.
parameter
A named vector: threshold
is the value that maximizes the Lagrange Multiplier values.
test.v
Vector of values of the two LM statistics for each threshold given in thd.range
.
thd.range
Range of values of the threshold.
fit.ARMA
The null model: ARMA fit over x
.
sigma2
Estimated innovation variance from the ARMA fit.
data.name
A character string giving the name of the data.
prop
Proportion of values of the series that fall in the lower regime.
p.value
The p-value of the test. It is NULL
for the asymptotic test.
method
A character string indicating the type of test performed.
d
The delay parameter.
pa
Lower threshold quantile.
dfree
Effective degrees of freedom. It is the number of tested parameters.
Simone Giannerini, [email protected]
Greta Goracci, [email protected]
Goracci G, Giannerini S, Chan K, Tong H (2023). “Testing for threshold effects in the TARMA framework.” Statistica Sinica, 33(3), 1879-1901. doi:10.5705/ss.202021.0120.
Andrews DWK (2003). “Tests for Parameter Instability and Structural Change with Unknown Change Point: A Corrigendum.” Econometrica, 71(1), 395-397. doi:10.1111/1468-0262.00405.
TAR.test
for the AR vs TAR asymptotic version of the test with different defaults. TAR.test.B
for the bootstrap version of the AR vs TAR test. TARMAGARCH.test
for the robust version of the test that assumes GARCH innovations. TARMA.sim
to simulate from a TARMA process.
## a TARMA(1,1,1,1) where the threshold effect is on the AR parameters set.seed(123) x1 <- TARMA.sim(n=100, phi1=c(0.5,-0.5), phi2=c(0.0,0.8), theta1=0.5, theta2=0.5, d=1, thd=0.2) TARMA.test(x1, ar.ord=1, ma.ord=1, d=1) TARMA.test(x1, ar.ord=1, ma.ord=1, d=1, ma.fixed=FALSE) # full TARMA test ## a TARMA(1,1,1,1) where the threshold effect is on the MA parameters set.seed(212) x2 <- TARMA.sim(n=100, phi1=c(0.5,0.2), phi2=c(0.5,0.2), theta1=0.6, theta2=-0.6, d=1, thd=0.2) TARMA.test(x2, ar.ord=1, ma.ord=1, d=1) TARMA.test(x2, ar.ord=1, ma.ord=1, d=1, ma.fixed=FALSE) # full TARMA test ## a ARMA(1,1) x3 <- arima.sim(n=100, model=list(order = c(1,0,1),ar=0.5, ma=0.5)) TARMA.test(x3, ar.ord=1, ma.ord=1, d=1) ## a TAR(1,1) x4 <- TARMA.sim(n=100, phi1=c(0.5,-0.5), phi2=c(0.0,0.8), theta1=0, theta2=0, d=1, thd=0.2) TARMA.test(x4, ar.ord=1, ma.ord=0, d=1) ## a AR(1) x5 <- arima.sim(n=100, model=list(order = c(1,0,0),ar=0.5)) TARMA.test(x5, ar.ord=1, ma.ord=0, d=1)
## a TARMA(1,1,1,1) where the threshold effect is on the AR parameters set.seed(123) x1 <- TARMA.sim(n=100, phi1=c(0.5,-0.5), phi2=c(0.0,0.8), theta1=0.5, theta2=0.5, d=1, thd=0.2) TARMA.test(x1, ar.ord=1, ma.ord=1, d=1) TARMA.test(x1, ar.ord=1, ma.ord=1, d=1, ma.fixed=FALSE) # full TARMA test ## a TARMA(1,1,1,1) where the threshold effect is on the MA parameters set.seed(212) x2 <- TARMA.sim(n=100, phi1=c(0.5,0.2), phi2=c(0.5,0.2), theta1=0.6, theta2=-0.6, d=1, thd=0.2) TARMA.test(x2, ar.ord=1, ma.ord=1, d=1) TARMA.test(x2, ar.ord=1, ma.ord=1, d=1, ma.fixed=FALSE) # full TARMA test ## a ARMA(1,1) x3 <- arima.sim(n=100, model=list(order = c(1,0,1),ar=0.5, ma=0.5)) TARMA.test(x3, ar.ord=1, ma.ord=1, d=1) ## a TAR(1,1) x4 <- TARMA.sim(n=100, phi1=c(0.5,-0.5), phi2=c(0.0,0.8), theta1=0, theta2=0, d=1, thd=0.2) TARMA.test(x4, ar.ord=1, ma.ord=0, d=1) ## a AR(1) x5 <- arima.sim(n=100, model=list(order = c(1,0,0),ar=0.5)) TARMA.test(x5, ar.ord=1, ma.ord=0, d=1)
Implements a supremum Lagrange Multiplier test for a ARMA-GARCH specification versus a TARMA-GARCH specification. Both the AR and MA parameters are tested. Also includes the ARCH case.
TARMAGARCH.test( x, pa = 0.25, pb = 0.75, ar.ord = 1, ma.ord = 1, arch.ord = 1, garch.ord = 1, d = 1, thd.range, ... )
TARMAGARCH.test( x, pa = 0.25, pb = 0.75, ar.ord = 1, ma.ord = 1, arch.ord = 1, garch.ord = 1, d = 1, thd.range, ... )
x |
A univariate time series. |
pa |
Real number in |
pb |
Real number in |
ar.ord |
Order of the AR part. It must be a positive integer |
ma.ord |
Order of the MA part. It must be a positive integer |
arch.ord |
Order of the ARCH part. It must be a positive integer |
garch.ord |
Order of the GARCH part. It must be a non negative integer. |
d |
Delay parameter. Defaults to |
thd.range |
Vector of optional user defined threshold range. If missing then |
... |
Additional arguments to be passed to |
Implements an asymptotic supremum Lagrange Multiplier test to test an ARMA-GARCH specification versus
a TARMA-GARCH specification. In other words, the test is robust with respect to heteroskedasticity.
Both the AR parameters and the MA parameters are tested. This is an asymptotic test and the value of
the test statistic has to be compared with the critical values reported in the output and taken from
(Andrews 2003). It includes the ARCH case if garch.ord=0
.
The null ARMA-GARCH model is estimated in one step with the function ugarchfit
from the package rugarch
. The estimated AR and MA polynomials are checked for stationarity and invertibility.
A list of class htest
with components:
statistic
The value of the supLM statistic.
parameter
A named vector: threshold
is the value that maximises the Lagrange Multiplier values.
test.v
Vector of values of the LM statistic for each threshold given in thd.range
.
thd.range
Range of values of the threshold.
fit
The null model: ARMA-GARCH fit using rugarch
.
sigma2
Estimated innovation variance from the ARMA fit.
data.name
A character string giving the name of the data.
prop
Proportion of values of the series that fall in the lower regime.
p.value
The p-value of the test. It is NULL
for the asymptotic test.
method
A character string indicating the type of test performed.
d
The delay parameter.
pa
Lower threshold quantile.
dfree
Effective degrees of freedom. It is the number of tested parameters.
Simone Giannerini, [email protected]
Greta Goracci, [email protected]
Angelini F, Castellani M, Giannerini S, Goracci G (2023). “Testing for Threshold Effects in Presence of Heteroskedasticity and Measurement Error with an application to Italian Strikes.” University of Bologna and Free University of Bolzano. 2308.00444, https://arxiv.org/abs/2308.00444.
Goracci G, Giannerini S, Chan K, Tong H (2021). “Testing for threshold effects in the TARMA framework.” University of Bologna, Free University of Bolzano, University of Iowa, London School of Economics. https://arxiv.org/abs/2103.13977.
Andrews DWK (2003). “Tests for Parameter Instability and Structural Change with Unknown Change Point: A Corrigendum.” Econometrica, 71(1), 395-397. doi:10.1111/1468-0262.00405.
TARMA.test
and TAR.test.B
for the asymptotic and bootstrap test without the GARCH component. TARMA.sim
to simulate from a TARMA process. TARMA.fit
and TARMA.fit2
for TARMA modelling.
## Function to simulate from a ARMA-GARCH process arma11.garch11 <- function(n, ph, th, a, b, a0=1, rand.gen = rnorm, innov = rand.gen(n, ...), n.start = 500, start.innov = rand.gen(n.start, ...),...){ # Simulates a ARMA(1,1)-GARCH(1,1) process # with parameters ph, th, a, b, a0. # x[t] <- ph*x[t-1] + th*eps[t-1] + eps[t] # eps[t] <- e[t]*sqrt(v[t]) # v[t] <- a0 + a*eps[t-1]^2 + b*v[t-1]; # ph : AR # th : MA # a : ARCH # b : GARCH # checks if(abs(a+b)>=1) stop("model is not stationary") if(b/(1-a)>=1) stop("model has infinite fourth moments") if (!missing(start.innov) && length(start.innov) < n.start) stop(gettextf("'start.innov' is too short: need %d points", n.start), domain = NA) e <- c(start.innov[1L:n.start], innov[1L:n]) ntot <- length(e) x <- v <- eps <- double(ntot) v[1] <- a0/(1.0-a-b); eps[1] <- e[1]*sqrt(v[1]) x[1] <- eps[1]; for(i in 2:ntot){ v[i] <- a0 + a*eps[i-1]^2 + b*v[i-1]; eps[i] <- e[i]*sqrt(v[i]) x[i] <- ph*x[i-1] + th*eps[i-1] + eps[i] } if (n.start > 0) x <- x[-(1L:n.start)] return(ts(x)); } ## ************************************************************************** ## Comparison between the robust and the non-robust test in presence of GARCH errors ## Simulates from the ARMA(1,1)-GARCH(1,1) set.seed(12) x1 <- arma11.garch11(n=100, ph=0.9, th=0.5, a=0.85, b=0.1, a0=1,n.start=500) TARMAGARCH.test(x1, ar.ord=1, ma.ord=1, arch.ord=1, garch.ord=1, d=1) TARMA.test(x1, ar.ord=1, ma.ord=1, d=1, ma.fixed=FALSE) ## a TARMA(1,1,1,1) where the threshold effect is on the AR parameters set.seed(123) x2 <- TARMA.sim(n=100, phi1=c(0.5,-0.5), phi2=c(0.0,0.8), theta1=0.5, theta2=0.5, d=1, thd=0.2) TARMAGARCH.test(x2, ar.ord=1, ma.ord=1, d=1)
## Function to simulate from a ARMA-GARCH process arma11.garch11 <- function(n, ph, th, a, b, a0=1, rand.gen = rnorm, innov = rand.gen(n, ...), n.start = 500, start.innov = rand.gen(n.start, ...),...){ # Simulates a ARMA(1,1)-GARCH(1,1) process # with parameters ph, th, a, b, a0. # x[t] <- ph*x[t-1] + th*eps[t-1] + eps[t] # eps[t] <- e[t]*sqrt(v[t]) # v[t] <- a0 + a*eps[t-1]^2 + b*v[t-1]; # ph : AR # th : MA # a : ARCH # b : GARCH # checks if(abs(a+b)>=1) stop("model is not stationary") if(b/(1-a)>=1) stop("model has infinite fourth moments") if (!missing(start.innov) && length(start.innov) < n.start) stop(gettextf("'start.innov' is too short: need %d points", n.start), domain = NA) e <- c(start.innov[1L:n.start], innov[1L:n]) ntot <- length(e) x <- v <- eps <- double(ntot) v[1] <- a0/(1.0-a-b); eps[1] <- e[1]*sqrt(v[1]) x[1] <- eps[1]; for(i in 2:ntot){ v[i] <- a0 + a*eps[i-1]^2 + b*v[i-1]; eps[i] <- e[i]*sqrt(v[i]) x[i] <- ph*x[i-1] + th*eps[i-1] + eps[i] } if (n.start > 0) x <- x[-(1L:n.start)] return(ts(x)); } ## ************************************************************************** ## Comparison between the robust and the non-robust test in presence of GARCH errors ## Simulates from the ARMA(1,1)-GARCH(1,1) set.seed(12) x1 <- arma11.garch11(n=100, ph=0.9, th=0.5, a=0.85, b=0.1, a0=1,n.start=500) TARMAGARCH.test(x1, ar.ord=1, ma.ord=1, arch.ord=1, garch.ord=1, d=1) TARMA.test(x1, ar.ord=1, ma.ord=1, d=1, ma.fixed=FALSE) ## a TARMA(1,1,1,1) where the threshold effect is on the AR parameters set.seed(123) x2 <- TARMA.sim(n=100, phi1=c(0.5,-0.5), phi2=c(0.0,0.8), theta1=0.5, theta2=0.5, d=1, thd=0.2) TARMAGARCH.test(x2, ar.ord=1, ma.ord=1, d=1)
Implements a supremum Lagrange Multiplier unit root test for the null hypiythesis of a integrated MA process versus a stationary TARMA process.
TARMAur.test(x, pa = 0.25, pb = 0.75, thd.range, method = "ML", ...)
TARMAur.test(x, pa = 0.25, pb = 0.75, thd.range, method = "ML", ...)
x |
A univariate vector or time series. |
pa |
Real number in |
pb |
Real number in |
thd.range |
Vector of optional user defined threshold range. If missing then |
method |
Fitting method to be passed to |
... |
Additional arguments to be passed to |
Implements an asymptotic supremum Lagrange Multiplier test to test an integrate MA(1) specification versus a
stationary TARMA(1,1) specification. This is an asymptotic test and the value of the test statistic has to be compared with the critical
values tabulated in (Chan et al. 2024) and available in supLMQur
.
The relevant critical values are automatically shown upon printing the test, see print.TARMAtest
.
An object of class TARMAtest
with components:
statistic
The value of the supLM statistic.
parameter
A named vector: threshold
is the value that maximises the Lagrange Multiplier values.
test.v
Vector of values of the LM statistic for each threshold given in thd.range
.
thd.range
Range of values of the threshold.
fit.ARMA
The null model: IMA(1) fit over x
.
sigma2
Estimated innovation variance from the IMA fit.
data.name
A character string giving the name of the data.
p.value
The p-value of the test. It is NULL
for the asymptotic test.
method
A character string indicating the type of test performed.
d
The delay parameter.
pa
Lower threshold quantile.
Simone Giannerini, [email protected]
Greta Goracci, [email protected]
Chan K, Giannerini S, Goracci G, Tong H (2024). “Testing for threshold regulation in presence of measurement error.” Statistica Sinica, 34(3), 1413-1434. doi:10.5705/ss.202022.0125, https://doi.org/10.5705/ss.202022.0125.
TARMAur.test.B
for the bootstrap version of the test.
print.TARMAtest
for the print method.
## a TARMA(1,1,1,1) set.seed(123) x1 <- TARMA.sim(n=100, phi1=c(0.5,-0.5), phi2=c(0.0,0.8), theta1=0.5, theta2=0.5, d=1, thd=0.2) TARMAur.test(x1) ## a IMA(1,1) x2 <- arima.sim(n=100, model=list(order = c(0,1,1),ma=0.6)) TARMAur.test(x2)
## a TARMA(1,1,1,1) set.seed(123) x1 <- TARMA.sim(n=100, phi1=c(0.5,-0.5), phi2=c(0.0,0.8), theta1=0.5, theta2=0.5, d=1, thd=0.2) TARMAur.test(x1) ## a IMA(1,1) x2 <- arima.sim(n=100, model=list(order = c(0,1,1),ma=0.6)) TARMAur.test(x2)
Implements a supremum Lagrange Multiplier unit root test for the null hypothesis of a integrated MA process versus a stationary TARMA process.
TARMAur.test.B( x, B = 1000, pa = 0.25, pb = 0.75, thd.range, method = "ML", btype = c("wb.r", "wb.n", "iid"), ... )
TARMAur.test.B( x, B = 1000, pa = 0.25, pb = 0.75, thd.range, method = "ML", btype = c("wb.r", "wb.n", "iid"), ... )
x |
A univariate vector or time series. |
B |
Integer. Number of bootstrap resamples. Defaults to 1000. |
pa |
Real number in |
pb |
Real number in |
thd.range |
Vector of optional user defined threshold range. If missing then |
method |
Fitting method to be passed to |
btype |
Bootstrap type, can be one of |
... |
Additional arguments to be passed to |
Implements the bootstrap version of TARMAur.test
the supremum Lagrange Multiplier test
to test an integrate MA(1) specification versus a stationary TARMA(1,1) specification.
The option btype
specifies the type of bootstrap as follows:
wb.r
Residual wild bootstrap with Rademacher auxiliary distribution. See (Giannerini et al. 2022).
wb.n
Residual wild bootstrap with Normal auxiliary distribution. See (Giannerini et al. 2022).
iid
Residual iid bootstrap. See (Goracci et al. 2021).
An object of class TARMAtest
with components:
statistic
The value of the supLM statistic.
parameter
A named vector: threshold
is the value that maximises the Lagrange Multiplier values.
test.v
Vector of values of the LM statistic for each threshold given in thd.range
.
thd.range
Range of values of the threshold.
fit.ARMA
The null model: IMA(1) fit over x
.
sigma2
Estimated innovation variance from the IMA fit.
data.name
A character string giving the name of the data.
p.value
The bootstrap p-value of the test.
method
A character string indicating the type of test performed.
d
The delay parameter.
pa
Lower threshold quantile.
Tb
The bootstrap null distribution.
Simone Giannerini, [email protected]
Greta Goracci, [email protected]
Chan K, Giannerini S, Goracci G, Tong H (2024). “Testing for threshold regulation in presence of measurement error.” Statistica Sinica, 34(3), 1413-1434. doi:10.5705/ss.202022.0125, https://doi.org/10.5705/ss.202022.0125.
TARMAur.test
for the asymptotic version of the test. print.TARMAtest
for the print method.
## a TARMA(1,1,1,1) set.seed(123) x1 <- TARMA.sim(n=100, phi1=c(0.5,-0.5), phi2=c(0.0,0.8), theta1=0.5, theta2=0.5, d=1, thd=0.2) TARMAur.test.B(x1, B=100) # B=100 for speedup ## a IMA(1,1) x2 <- arima.sim(n=100, model=list(order = c(0,1,1),ma=0.6)) TARMAur.test.B(x2, B=100) # B=100 for speedup
## a TARMA(1,1,1,1) set.seed(123) x1 <- TARMA.sim(n=100, phi1=c(0.5,-0.5), phi2=c(0.0,0.8), theta1=0.5, theta2=0.5, d=1, thd=0.2) TARMAur.test.B(x1, B=100) # B=100 for speedup ## a IMA(1,1) x2 <- arima.sim(n=100, model=list(order = c(0,1,1),ma=0.6)) TARMAur.test.B(x2, B=100) # B=100 for speedup