Package 'tsDyn'

Description

Getting started with the tsDyn package

Details

This package provide some tools inspired by nonlinear dynamics for the analysis-modelling of observed time series.
For loading the package, type:
library(tsDyn)
A good place to start learning the package usage, is the vignette. It contains a more detailed guide on package contents, and an applied case study. At the R prompt, write:
vignette("tsDyn")
For a full list of functions exported by the package, type:
ls("package:tsDyn")

There is also an experimental GUI for built-in NLAR models. Call it with:
nlarDialog(timeSeries)
where timeSeries is an available time series object.

Each exported function has a corresponding man page (some man pages are in common to more functions). Display it by typing
help(functionName)

Author(s)

Antonio, Fabio Di Narzo

See Also

availableModels for listing all currently available NLAR models

autopairs,autotriples,autotriples.rgl for graphical exploratory functions

llar, delta, delta.lin for nonlinearity checking tools

Description

Additive nonlinear autoregressive model.

Usage

aar(x, m, d=1, steps=d, series)

Arguments

Details

Nonparametric additive autoregressive model of the form:

$x_{t+s} = \mu + \sum_{j=1}^{m} s_j(x_{t-(j-1)d})$

where $s_j$ are nonparametric univariate functions of lagged time series values. They are represented by cubic regression splines. $s_j$ are estimated together with their level of smoothing using routines in the mgcv package (see references).

Value

An object of class nlar, subclass aar, i.e. a list with mostly internal structures for the fitted gam object.

Author(s)

Antonio, Fabio Di Narzo

References

Wood, mgcv:GAMs and Generalized Ridge Regression for R. R News 1(2):20-25 (2001)

Wood and Augustin, GAMs with integrated model selection using penalized regression splines and applications to environmental modelling. Ecological Modelling 157:157-177 (2002)

Examples

#fit an AAR model:
mod <- aar(log(lynx), m=3)
#Summary informations:
summary(mod)
#Diagnostic plots:
plot(mod)

Description

Compute forecasting accuracies. This is very similar to the accuracy method form forecast.

Usage

accuracy_stat(object, ...)

## Default S3 method:
accuracy_stat(object, true, ...)

## S3 method for class 'pred_roll'
accuracy_stat(object, w, ...)

Arguments

Details

The function works either for a simple data.frame or for objects pred_roll. For simple data.frames, the argument true, i.e. a data frame containing the true values, has to be provided. For pred_roll objects, the true values are contained in the object, so no need (nor possibility) to provide the true values.

Value

A data-frame containing the forecasting accuracy measures.

Author(s)

Matthieu Stigler

Examples

## univariate:
mod_ar <- linear(lynx[1:100], m=1)
mod_ar_pred <- predict_rolling(mod_ar, newdata=lynx[101:114])
accuracy_stat(object=mod_ar_pred$pred, true=mod_ar_pred$true)

## multivariate
data(barry)
mod_var <- lineVar(barry, lag=1)

mod_var_pred <-predict_rolling(object=mod_var, nroll=10, n.ahead=1:3)
accuracy_stat(object=mod_var_pred)
accuracy_stat(object=mod_var_pred, w=c(0.7, 0.2, 0.1))

Description

addRegime test

Usage

addRegime(object, ...)

Arguments

Value

A list containing the p-value of the F statistic and a boolean, true if there is some remaining nonlinearity and false otherwise.

Author(s)

J. L. Aznarte

References

TODO

See Also

star

Examples

##TODO

Description

Computes the long term mean of an AR process

Usage

ar_mean(object, ...)

## S3 method for class 'linear'
ar_mean(object, ...)

## S3 method for class 'setar'
ar_mean(object, ...)

## S3 method for class 'lstar'
ar_mean(object, ...)

Arguments

Details

The function computes the long-term mean of an AR(p) process, or of the corresponding sub-regimes in SETAR or LSTAR model. There are three possible cases:

No constant nor trend

The LT mean is 0

constant

The LT mean is given by const/(1-sum(AR coefs))

Trend

The LT mean is not defined

Examples

## estimate a (linear) AR, a SETAR and a LSTAR
lin_cst_l1 <-  linear(lh, m = 1, include = "const")
set_cst_l1 <-  setar(lh, m = 1, include = "const") 
lst_cst_l1 <-  lstar(lh, m = 1, include = "const", trace = FALSE)

ar_mean(lin_cst_l1)
ar_mean(set_cst_l1)
ar_mean(lst_cst_l1)

Description

Bivariate time series plots: scatterplots, directed lines and kernel density estimations using functions in the sm package.

Usage

autopairs(
  x,
  lag = 1,
  h,
  type = c("levels", "persp", "image", "lines", "points", "regression")
)

Arguments

Details

Bivariate time series plots: scatterplots, directed lines and kernel density and regression functions estimations using functions in the package sm. In particular, for kernel density estimation sm.density is used, with smoothing parameter h defaulting to hnorm. For kernel regression, sm.regression is used.

Value

None. Plots are produced on the default graphical device.

Author(s)

Wrappers to sm by Antonio, Fabio Di Narzo

See Also

For finer control on density estimation, consider using directly sm.density and, especially, sm.ts.pdf from package sm.

Examples

x <- log10(lynx)
autopairs(x, lag=2, type="lines")

Description

Trivariate time series plots: kernel autoregression using functions in the sm package

Usage

autotriples(
  x,
  lags = 1:2,
  h,
  type = c("levels", "persp", "image", "lines", "points")
)

Arguments

Details

This function displays trivariate time series plots, i.e. kernel regression of $x[t-lags[1]], x[t-lags[2]]$ against $x[t]$ using functions in the package sm. In particular, sm.regression is used, with smoothing parameter defaulting to hnorm(x).

Value

None. Plots are produced on the default graphical device.

Author(s)

Wrappers to sm by Antonio, Fabio Di Narzo

See Also

For finer control on kernel regression, consider using directly sm.regression and, especially, sm.autoregression in package sm.

Examples

autotriples(log(lynx))
autotriples(log(lynx), type="persp")
autotriples(log(lynx), type="image")

Description

Interactive trivariate time series plots

Usage

autotriples.rgl(x, lags = 1:2, type = c("lines", "points"))

Arguments

Details

This function displays interactive trivariate time series plots x[t-lags[1]], x[t-lags[2]] against x[t] using the interactive rgl device.

Value

None. A plot is produced on the current rgl device.

Author(s)

Wrapper to 'sm' and GUI by Antonio, Fabio Di Narzo

See Also

autotriples for 3d visualization via scatterplot3d package and for kernel post-processing of the cloud for nonparametric autoregression functions estimates.

Examples

if(interactive())
autotriples.rgl(log(lynx))

Description

Available built-in time series models

Usage

availableModels()

Details

Return the list of built-in available ‘nlar’ time series models

Value

A character vector containing built-in time series models. For help on a specific model, type: help(modelName).

Author(s)

Antonio, Fabio Di Narzo

Examples

availableModels()

Description

This data set contains the series used by Bierens and Martins for testing for PPI between Canada and US.

Usage

data(barry)

Format

A data frame with 324 monthly observations, ranging from 1973:M1 until 1999:M12.

Author(s)

Matthieu Stigler

Source

Bierens, H. and Martins, L. (2010), Time Varying Cointegration,

Description

Test of unit root against a stationnary three regime SETAR alternative

Usage

BBCTest(
  x,
  m,
  series,
  testStat = c("LR", "Wald", "LM"),
  trim = 0.1,
  grid = c("minPerc", "minObs")
)

Arguments

Details

TODO

Value

A object of class "BBC2004Test" containing:

-The value of the sup Test

-The version of test used (either Wald, LM or LR).

Author(s)

Matthieu Stigler

References

Bec, Ben Salem and Carrasco (2004) Tests for Unit-Root versus Threshold Specification With an Application to the Purchasing Power Parity Relationship, Journal of Business and Economic Statistics; 22:4

See Also

setarTest for a test with stationarity as a null.

Examples

BBCTest(lynx, m=3, test="Wald", grid="minPerc")

Description

Computes the (inverse) characteristic roots of the auto-regressive coefficients. To be stationary, the values should be outside the unit circle. The function here returns the modulus of the roots.

Usage

charac_root(object, ...)

## S3 method for class 'nlar'
charac_root(object, ...)

Arguments

Details

Computes the roots of the polynomial (1, -phi) using function polyroot

Value

a data.frame, with the modulus of the roots. For models with multiple regimes (setar, lstar, star), one column per regime.

Examples

mod.ar <-  linear(lh, m = 5, include = "const")
mod.setar <-  setar(lh, m = 5, include = "const")

charac_root(mod.ar)
charac_root(mod.setar)

Description

Extract parameters in VECM: adjustment coefficients A, cointegrating coefficients B , or the composite matrix PI

Usage

coefB(object, ...)

## S3 method for class 'VECM'
coefB(object, ...)

## S3 method for class 'ca.jo'
coefB(object, r = 1, normalize = TRUE, ...)

coefA(object, ...)

## S3 method for class 'VECM'
coefA(object, ...)

## S3 method for class 'ca.jo'
coefA(object, r = 1, normalize = TRUE, ...)

coefPI(object, ...)

Arguments

Details

The functions extract the parameters from a VECM with $K$ variables and rank $r$:

A

Adjustment coefficients, of dim $K \times r$

B

Cointegrating coefficients, of dim $K \times r$

Pi

Matrix $\Pi=A B^{'}$, of dim $K \times K$

Coefficients are extracted from a VECM in package tsDyn, or from a VECM obtained in package urca from ca.jo or cajorls.

Note that by default, the A and B coefficients returned are normalized (see below). This is the case for results obtained from VECM/lineVar and cajorls, while for ca.jo, the user has the choice (but normalize=TRUE by default), in which case the rank r is also to be specified. The normalization is the Phillips triangular representation, as suggested by Johansen (1995, p. 72), standardising the first $r\times r$ coefficients to $I_r$:

B

$B_{norm}=B (c^{'}B)^{-1}$ with $c=(I_r,0_{K-r,r})^{'}$

A

$A_{norm}=AB^{'}c$

Finally, note that the function also apply to objects obtained from tests of class ca.jo.test (from blrtest etc...). Care should be taken however, since the normalization might override the restrictions imposed.

Value

A matrix containing the coefficients

Author(s)

Matthieu Stigler

References

Johansen, Soren, (1995), Likelihood-Based Inference in Cointegrated Vector Autoregressive Models, Oxford University Press

Examples

data(barry)
vecm <- VECM(barry,  lag=1, estim="ML")
vecm_r2 <- VECM(barry,  lag=1, estim="ML", r=2)

## extract coefficients:
coefA(vecm)
coefB(vecm)
coefPI(vecm)
coefB(vecm_r2)
coefPI(vecm_r2)

## Beta-Restricted VECM:
beta_vecm2 <- coefB(vecm_r2) 
beta_vecm2[3,2] <- 0.02
vecm_r2_rest <- VECM(barry,  lag=1, estim="ML", r=2, beta=beta_vecm2)
round(coefB(vecm_r2_rest),5)

## Package vars/urca
if(require(urca)){
 vecm_ur <- ca.jo(barry, K=2)
 coefB(vecm_ur)
 coefB(vecm_ur,r=2)
 coefB(cajorls(vecm_ur, r=2))
 all.equal(coefB(vecm), coefB(vecm_ur), check.attributes=FALSE)
 all.equal(coefB(vecm_r2), coefB(vecm_ur, r=2), check.attributes=FALSE)
}

Description

delta statistic of conditional independence and associated bootstrap test

Usage

delta(x, m, d = 1, eps)

delta.test(
  x,
  m = 2:3,
  d = 1,
  eps = seq(0.5 * sd(x), 2 * sd(x), length = 4),
  B = 49
)

Arguments

Details

delta statistic of conditional independence and associated bootstrap test. For details, see Manzan(2003).

Value

delta returns the computed delta statistic. delta.test returns the bootstrap based 1-sided p-value.

Warning

Results are sensible to the choice of the window eps. So, try the test for a grid of m and eps values. Also, be aware of the course of dimensionality: m can't be too high for relatively small time series. See references for further details.

Author(s)

Antonio, Fabio Di Narzo

References

Sebastiano Manzan, Essays in Nonlinear Economic Dynamics, Thela Thesis (2003)

See Also

BDS marginal independence test: bds.test in package tseries

Teraesvirta's neural network test for nonlinearity: terasvirta.test in package tseries

delta test for nonlinearity: delta.lin.test

Examples

delta(log10(lynx), m=3, eps=sd(log10(lynx)))

Description

delta test of linearity based on conditional mutual information

Usage

delta.lin(x, m, d = 1)

delta.lin.test(
  x,
  m = 2:3,
  d = 1,
  eps = seq(0.5 * sd(x), 2 * sd(x), length = 4),
  B = 49
)

Arguments

Details

delta test of linearity based on conditional mutual information

Value

delta.lin returns the parametrically estimated delta statistic for the given time series (assuming linearity). delta.lin.test returns the bootstrap based 1-sided p-value. The test statistic is the difference between the parametric and nonparametric delta estimators.

Author(s)

Antonio, Fabio Di Narzo

References

Sebastiano Manzan, Essays in Nonlinear Economic Dynamics, Thela Thesis (2003)

Examples

delta.lin(log10(lynx), m=3)

Description

Use the fevd function from package vars to compute the forecast error variance decomposition of a VAR(p) or VECM for n.ahead steps.

Usage

## S3 method for class 'nlVar'
fevd(x, n.ahead = 10, ...)

Arguments

Details

The function converts the VAR or VECM computed by package tsDyn into an object of class ‘vec2var’, on which then the fevd method is applied. For details, see package vars.

Value

A list with class attribute ‘varfevd’ of length K holding the forecast error variances as matrices.

Author(s)

Bernhard Pfaff

References

Hamilton, J. (1994), Time Series Analysis, Princeton University Press, Princeton.

Lutkepohl, H. (2006), New Introduction to Multiple Time Series Analysis, Springer, New York.

See Also

plot for the plot method. lineVar, VECM for the models.

Examples

data(zeroyld)
mod_vecm <- VECM(zeroyld, lag = 2)
fevd(mod_vecm, n.ahead = 5)

Description

Returns the fitted values of the model, either as computed in the model, or back to the original series level.

Usage

## S3 method for class 'nlVar'
fitted(object, level = c("model", "original"), ...)

Arguments

Details

In case of a VAR in differences, in ADF specification, or a VECM, the fitted values are actually in differences. With the option level="original", the function returns the series in the original level.

For VAR in levels, the two arguments are evidently the same and hence it is not taken into account, returning a warning.

Value

A matrix.

Author(s)

Matthieu Stigler

Examples

## estimate models
data(barry)

ve <- VECM(barry, lag=2)
va <- lineVar(barry, lag=1)
va_diff <- lineVar(barry, lag=1, I="diff")
va_ADF <- lineVar(barry, lag=1, I="ADF")


## get fitted values:
tail(fitted(ve))
tail(fitted(ve, level="original"))

tail(fitted(va))
tail(fitted(object=va, level="original"))

tail(fitted(va_diff))
tail(fitted(object=va_diff, level="original"))

tail(fitted(va_ADF))
tail(fitted(object=va_ADF, level="original"))

Description

Extract threshold coefficient(s)

Usage

getTh(object, ...)

## Default S3 method:
getTh(object, ...)

Arguments

Value

Threshold value.

Author(s)

Matthieu Stigler

Examples

set<-setar(lynx, m=3)
getTh(set)
getTh(summary(set))

Description

Generates a GIRF for multiple innovations and histories

Usage

GIRF(object, n.ahead, seed = NULL, ...)

## S3 method for class 'setar'
GIRF(
  object,
  n.ahead = 10,
  seed = NULL,
  n.hist = 20,
  n.shock = 20,
  R = 10,
  hist_li = NULL,
  shock_li = NULL,
  ...
)

## S3 method for class 'linear'
GIRF(
  object,
  n.ahead = 10,
  seed = NULL,
  n.hist = 20,
  n.shock = 20,
  R = 10,
  hist_li = NULL,
  shock_li = NULL,
  ...
)

## S3 method for class 'nlVar'
GIRF(
  object,
  n.ahead = 10,
  seed = NULL,
  n.hist = 20,
  n.shock = 20,
  R = 10,
  hist_li = NULL,
  shock_li = NULL,
  ...
)

## S3 method for class 'GIRF_df'
plot(
  x,
  plot_type = c("density", "line"),
  n.ahead = c(1, 5, 10),
  var = unique(x$var)[1],
  n_simu = c(1, 2),
  ...
)

Arguments

Details

In a nonlinear model, the Impulse response Function (IRF) is not time-, scale- and sign-invariant as in linear models. To address this, Koop et al (1996) introduced the Generalized Impulse response Function (GIRF):

$GIRF(h,\delta,\omega_{t-1})=E[y_{t+h}|\epsilon_{t}=\delta,\epsilon_{t+h}=0,\omega_{t-1}]-E[y_{t+h}|\epsilon_{t}=0,\epsilon_{t+h}=0,\omega_{t-1}]$

It is the difference between two conditional expectations, one containing the shock of interest, the second one averaging it out. The averaging-out is done by comparing against random innovations (unlike the IRF, which compares against innovation 0), The parameter R corresponds to the number of times this is done.

The GIRF as defined here depends on the particular shock, as well as on the history. Koop et al (1996) suggest to draw multiple combinations of histories and innovations. This is done with arguments n.hist and n.shock (or, alternatively, provide one of, or both, hist_li and hist_li as list of histories and shocks).

The output is a data-frame containing the two average paths for each combinations of shocks and histories.

Value

A data-frame, with:

n_simu:

Id for the simulation (total number is n.hist times n.shock)

hist, shock

History and shock used in the nth simulation

n.ahead:

The forecasting horizon. Note the shocks happens at time 0

var:

The variable (on which the shock happens, corresponds hence to the response argument in irf)

sim_1, sim_2

The average (over R times) simulation with the specific shock (sim_1) or with random shocks (sim_2).

girf

The difference between sim_1 and sim_2

Author(s)

Matthieu Stigler

References

Koop, G, Pesaran, M. H. & Potter, S. M. (1996) Impulse response analysis in nonlinear multivariate models. Journal of Econometrics, 74, 119-147

See Also

irf.nlVar for the IRF, for linear models, or in case of non-linear models, for each regime.

Examples

## simulate a SETAR for the example. Higher regime more persistent (AR coef 0.5 instead of 0.2)
set <- setar.sim(B = c(0.5, 0.2, 0.5, 0.5), lag = 1, Thresh = 0.5, n = 500)
set_estim <- setar(set, m = 1)

## regime-specific IRF:
plot(irf(set_estim, regime = "L", boot = FALSE))
plot(irf(set_estim, regime = "H", boot = FALSE))

## GIRF
girf_out <- GIRF(set_estim, n.hist = 10, n.shock = 10) # smaller number for example only

## the GIRF shows a very fast convergence (the shock at n.ahead = 4 is already very close to 0)
plot(girf_out, n.ahead = 1:4)
## investigate a few specific GIRFS:
plot(girf_out, plot_type = "line", n_simu  = 1:5)

Description

This data, used as example in Hansen (1999), contains the US monthly industrial production.

Usage

IIPUs

Format

A monthly time series of class ts starting in January 1960 and ending in September 1997. Note that the series ends at 1997 and not 1998 as in the paper of Hansen, even if the data was taken from hi site and the graph is exactly the same.

Source

Hansen (1999) Testing for linearity, Journal of Economic Surveys, Volume 13, Number 5, December 1999 , pp. 551-576(26) available at: https://www.ssc.wisc.edu/~bhansen/papers/cv.htm

Examples

data(IIPUs)
end(IIPUs) #not same date as in the paper

## Figure 2 in paper (p. 555)
plot(IIPUs)

## Table 1 in paper (p. 556)
ar_1 <- linear(IIPUs, m=16)
coef(summary(ar_1))[, 1:2]
deviance(ar_1)

## Table 2 in paper (p. 559)
set_1 <- setar(IIPUs, m=16, thDelay=5, trim=0.1)
## tsDyn finds another threshold, with a better SSR:
getTh(set_1)
deviance(set_1)

## estimate with Hansen threshold:
set_1_han <- setar(IIPUs, m=16, thDelay=5, trim=0.1, th = 0.23)
deviance(set_1_han)
set_1_CO <- coef(summary(set_1_han))[, 1:2]
cbind(set_1_CO[1:17,], set_1_CO[18:34,])


## Table 4 in paper (p. 561)
set_2 <- setar(IIPUs, m=16, thDelay=5, trim=0.1, nthresh =2)
getTh(set_2)
deviance(set_2)

##most of the results agree, except constant in the low regime which has opposed signs!
set_2_CO <- coef(summary(set_2))[, 1:2]
cbind(set_2_CO[1:17,], set_1_CO[18:34,])

#this is obviously a error in Hansen, see:
XX<-embed(IIPUs, 17)
Y<-XX[,1]
X<-XX[,-1]
dummyDown<-ifelse(X[,6]<= -2.5, 1,0)
sum(dummyDown)
M<-cbind(1*dummyDown,X*dummyDown )
lm(Y~M-1)

## setar test (takes long to compute, even with small nboot)
## Not run: 
  test_1 <- setarTest(IIPUs, m=16, thDelay=5, nboot=10)
  #because of the discrepency. test1vs2 does not correspond, test 1vs3 conforms
  test_1$Ftests ## compare with Table 2 (F12)  and Table 4 (F13, F23)
  summary(test_1)
  
  setarTest(IIPUs, m=16, thDelay=5, nboot=10, test="2vs3")
  #test 2vs3 is also different of the version in the article (27)

## End(Not run)

## results from the test is stored in: setarTest_IIPUs_results
data(setarTest_IIPUs_results)

## Table 5 and 6
test_1vs <- setarTest_IIPUs_results$test_1
test_1vs

## Table 7
test_2vs <- setarTest_IIPUs_results$test_2
test_2vs
plot(test_2vs)

Description

Compute the impulse response coefficients (IRF) of a VAR(p) (or transformed VECM to VAR(p)) for n.ahead steps. For TVECM and TVAR model, argument regime offers regime-specific IRF.

Usage

## S3 method for class 'linear'
irf(
  x,
  impulse = NULL,
  response = NULL,
  n.ahead = 10,
  ortho = TRUE,
  cumulative = FALSE,
  boot = TRUE,
  ci = 0.95,
  runs = 100,
  seed = NULL,
  ...
)

## S3 method for class 'setar'
irf(
  x,
  impulse = NULL,
  response = NULL,
  n.ahead = 10,
  ortho = TRUE,
  cumulative = FALSE,
  boot = TRUE,
  ci = 0.95,
  runs = 100,
  seed = NULL,
  regime = c("L", "M", "H"),
  ...
)

## S3 method for class 'ar'
irf(
  x,
  impulse = NULL,
  response = NULL,
  n.ahead = 10,
  ortho = TRUE,
  cumulative = FALSE,
  boot = TRUE,
  ci = 0.95,
  runs = 100,
  seed = NULL,
  ...
)

## S3 method for class 'VAR'
irf(
  x,
  impulse = NULL,
  response = NULL,
  n.ahead = 10,
  ortho = TRUE,
  cumulative = FALSE,
  boot = TRUE,
  ci = 0.95,
  runs = 100,
  seed = NULL,
  ...
)

## S3 method for class 'VECM'
irf(
  x,
  impulse = NULL,
  response = NULL,
  n.ahead = 10,
  ortho = TRUE,
  cumulative = FALSE,
  boot = TRUE,
  ci = 0.95,
  runs = 100,
  seed = NULL,
  ...
)

## S3 method for class 'TVAR'
irf(
  x,
  impulse = NULL,
  response = NULL,
  n.ahead = 10,
  ortho = TRUE,
  cumulative = FALSE,
  boot = TRUE,
  ci = 0.95,
  runs = 100,
  seed = NULL,
  regime = c("L", "M", "H"),
  ...
)

## S3 method for class 'TVECM'
irf(
  x,
  impulse = NULL,
  response = NULL,
  n.ahead = 10,
  ortho = TRUE,
  cumulative = FALSE,
  boot = TRUE,
  ci = 0.95,
  runs = 100,
  seed = NULL,
  regime = c("L", "M", "H"),
  ...
)

Arguments

Value

A list of class ‘varirf’ see irf.varest

Author(s)

Matthieu Stigler

See Also

plot for the plot method. lineVar, VECM for the models. Hamilton, J. (1994), Time Series Analysis, Princeton University Press, Princeton.

Lutkepohl, H. (2006), New Introduction to Multiple Time Series Analysis, Springer, New York.

Examples

data(barry)

## For VAR
mod_var <- lineVar(barry, lag = 2)
irf(mod_var, impulse = "dolcan", response = c("dolcan", "cpiUSA", "cpiCAN"), boot =
FALSE)

## For VECM
mod_VECM <- VECM(barry, lag = 2, estim="ML", r=2)
irf(mod_VECM, impulse = "dolcan", response = c("dolcan", "cpiUSA", "cpiCAN"), boot =
FALSE)

Description

Generic NLAR linearity test

Usage

isLinear(object, ...)

Arguments

Author(s)

A. F. Di Narzo

Description

Test of unit root against a stationnary 3 regime SETAR alternative with random walk in the inner regime

Usage

KapShinTest(x, m=1, series, include = c("none","const", "trend", "both"), 
  c=3, delta=0.5, points=NULL,minObsMid=10, 
  trick=c("for", "apply", "mapply"), trace=FALSE)

Arguments

Details

This function is currently spurious.

Value

A object of class KapShin2006Test containing:

Author(s)

Matthieu Stigler

See Also

BBCTest for a similar test. setarTest for a test with stationarity as a null.

Examples

KapShinTest(lynx, m=1, trace=FALSE, include="none", points=10)

Description

Selection of the cointegrating rank and the lags with Information criterion (AIC, BIC).

Usage

lags.select(
  data,
  lag.max = 10,
  include = c("const", "trend", "none", "both"),
  fitMeasure = c("SSR", "LL"),
  sameSample = TRUE
)

Arguments

Details

This function selects the lag according to AIC, BIC and Hannan-Quinn.

Value

An object of class rank.select, with ‘print’ and ‘summary methods’, containing among other the matrices of AIC/BIC/HQ.

Author(s)

Matthieu Stigler

See Also

rank.select, the underlying function, to estimate the rank also.

VARselect in package vars, does basically the same.

Examples

data(barry)

# 
rk_sel <- lags.select(barry)
rk_sel
summary(rk_sel)

Description

AR(m) model

Usage

linear(x, m, d=1, steps=d, series, include = c( "const", "trend","none", "both"), 
        type=c("level", "diff", "ADF"), warn_root=TRUE)

Arguments

Details

AR(m) model:

$x_{t+s} = \phi_0 + \phi_1 x_t + \phi_2 x_{t-d} + \dots + \phi_m x_{t - (m-1)d} + \epsilon_{t+s}$

Value

A nlar object, linear subclass.

Author(s)

Antonio, Fabio Di Narzo

See Also

nlar for fitting this and other models to time series data

Examples

#fit an AR(2) model
mod.linear <- linear(log(lynx), m=2)
mod.linear
summary(mod.linear)

Description

Estimate either a VAR or a VECM.

Usage

lineVar(
  data,
  lag,
  r = 1,
  include = c("const", "trend", "none", "both"),
  model = c("VAR", "VECM"),
  I = c("level", "diff", "ADF"),
  beta = NULL,
  estim = c("2OLS", "ML"),
  LRinclude = c("none", "const", "trend", "both"),
  exogen = NULL
)

Arguments

Details

This function provides basic functionalities for VAR and VECM models. More comprehensive functions are in package vars. A few differences appear in the VECM estimation:

Engle-Granger estimator

The Engle-Granger estimator is available

Presentation

Results are printed in a different ways, using a matrix form

lateX export

The matrix of coefficients can be exported to latex, with or without standard-values and significance stars

Two estimators are available: the Engle-Granger two-steps approach (2OLS) or the Johansen (ML). For the 2OLS, deterministic regressors (or external variables if LRinclude is of class numeric) can be added for the estimation of the cointegrating value and for the ECT. This is only working when the beta value is not pre-specified.

The argument beta is only for VECM, look at the specific help page for more details.

Value

Fitted model data

Author(s)

Matthieu Stigler

See Also

VECM which is just a wrapper for lineVar(..., model="VECM"). Methods predict.VAR, VARrep, regime, irf and toLatex.

TVAR and TVECM for the corresponding threshold models. linear for the univariate AR model.

Examples

data(zeroyld)

#Fit a VAR
VAR <- lineVar(zeroyld, lag=1)
VAR
summary(VAR)

#compare results with package vars:
if(require(vars)) {
a<-VAR(zeroyld, p=1)
coef_vars <- t(sapply(coef(a), function(x) x[c(3,1,2),1]))
all.equal(coef(VAR),coef_vars, check.attributes=FALSE)
}

###VECM
VECM.EG <- lineVar(zeroyld, lag=2, model="VECM")
VECM.EG
summary(VECM.EG)

VECM.ML <- lineVar(zeroyld, lag=2, model="VECM", estim="ML")
VECM.ML
summary(VECM.ML)


###Check Johansen MLE
myVECM <- lineVar(zeroyld, lag=1, include="const", model="VECM", estim="ML")
summary(myVECM, digits=7) 
#comparing with vars package
if(require(vars)){
a<-ca.jo(zeroyld, spec="trans")
summary(a)
#same answer also!
}

##export to Latex
toLatex(VECM.EG)
toLatex(summary(VECM.EG))
options("show.signif.stars"=FALSE)
toLatex(summary(VECM.EG), parenthese="Pvalue")
options("show.signif.stars"=TRUE)

Description

Casdagli test of nonlinearity via locally linear forecasts

Usage

llar(x, m, d = 1, steps = d, series, eps.min = sd(x)/2,
	eps.max = diff(range(x)), neps = 30, trace = 0)

llar.predict(x, m, d=1, steps=d, series, n.ahead=1, 
eps=stop("you must specify a window value"),
onvoid=c("fail","enlarge"), r = 20, trace=1)

llar.fitted(x, m, d=1, steps=d, series, eps, trace=0)

Arguments

Details

llar does the Casdagli test of non-linearity. Given the embedding state-space (of dimension m and time delay d) obtained from time series series, for a sequence of distance values eps, the relative error made by forecasting time series values with a linear autoregressive model estimated on points closer than eps is computed. If minimum error is reached at relatively small length scales, a global linear model may be inappropriate (using current embedding parameters). This was suggested by Casdagli(1991) as a test for non-linearity.

llar.predict tries to extend the given time series by n.ahead points by iteratively fitting locally (in the embedding space of dimension m and time delay d) a linear model. If the spatial neighbourhood window is too small, your time series last point would be probably isolated. You can ask to automatically enlarge the window eps by a factor of r% sequentially, until enough neighbours are found for fitting the linear model.

llar.fitted gives out-of-sample fitted values from locally linear models.

Value

llar gives an object of class 'llar'. I.e., a list of components:

Function llar.forecast gives the vector of n steps ahead locally linear iterated forecasts.

Function llar.fitted gives out-of-sample fitted values from locally linear models.

Warning

For long time series, this can be slow, especially for relatively big neighbourhood sizes.

Note

The C implementation was re-adapted from that in the TISEAN package ("ll-ar" routine, see references). However, here the euclidean norm is used, in place of the max-norm.

Author(s)

Antonio, Fabio Di Narzo

References

M. Casdagli, Chaos and deterministic versus stochastic nonlinear modelling, J. Roy. Stat. Soc. 54, 303 (1991)

Hegger, R., Kantz, H., Schreiber, T., Practical implementation of nonlinear time series methods: The TISEAN package; CHAOS 9, 413-435 (1999)

Examples

res <- llar(log(lynx), m=3, neps=7)
plot(res)

x.new <- llar.predict(log(lynx),n.ahead=20, m=3, eps=1, onvoid="enlarge", r=5)
lag.plot(x.new, labels=FALSE)

x.fitted <-  llar.fitted(log(lynx), m=3, eps=1)
lag.plot(x.fitted, labels=FALSE)

Description

Log-Likelihood method for VAR and VECM models.

Usage

## S3 method for class 'nlVar'
logLik(object, ...)

## S3 method for class 'VECM'
logLik(object, r, ...)

Arguments

Details

For a VAR, the Log-Likelihood is computed as in Luetkepohl (2006) equ. 3.4.5 (p. 89) and Juselius (2006) p. 56:

$LL = -(TK/2) \log(2\pi) - (T/2) \log|\Sigma| - (1/2) \sum^{T} \left [ (y_t - A^{'}x_t)^{'} \Sigma^{-1} (y_t - A^{'}x_t) \right ]$

Where $\Sigma$ is the Variance matrix of residuals, and $x_t$ is the matrix stacking the regressors (lags and deterministic).

However, we use a computationally simpler version:

$LL = -(TK/2) \log(2\pi) - (T/2) \log|\Sigma| - (TK/2)$

See Juselius (2006), p. 57.

(Note that Hamilton (1994) 11.1.10, p. 293 gives $+ (T/2) \log|\Sigma^{-1}|$, which is the same as $-(T/2) \log|\Sigma|)$.

For VECM, the Log-Likelihood is computed in two different ways, depending on whether the VECM was estimated with ML (Johansen) or 2OLS (Engle and Granger).

When the model is estimated with ML, the LL is computed as in Hamilton (1994) 20.2.10 (p. 637):

$LL = -(TK/2) \log(2\pi) - (TK/2) -(T/2) \log|\hat\Sigma_{UU}| - (T/2) \sum_{i=1}^{r} \log (1-\hat\lambda_{i})$

Where $\Sigma_{UU}$ is the variance matrix of residuals from the first auxiliary regression, i.e. regressing $\Delta y_t$ on a constant and lags, $\Delta y_{t-1}, \ldots, \Delta y_{t-p}$. $\lambda_{i}$ are the eigenvalues from the $\Sigma_{VV}^{-1}\Sigma_{VU}\Sigma_{UU}^{-1}\Sigma_{UV}$, see 20.2.9 in Hamilton (1994).

When the model is estimated with 2OLS, the LL is computed as:

$LL = \log|\Sigma|$

Where $\Sigma$ is the variance matrix of residuals from the the VECM model. There is hence no correspondence between the LL from the VECM computed with 2OLS or ML.

Value

Log-Likelihood value.

Author(s)

Matthieu Stigler

References

Hamilton (1994) Time Series Analysis, Princeton University Press

Juselius (2006) The Cointegrated VAR model: methodology and Applications, Oxford Univesity Press

Luetkepohl (2006) New Introduction to Multiple Time Series Analysis, Springer

Examples

data(zeroyld)

#Fit a VAR
VAR <- lineVar(zeroyld, lag=1)
logLik(VAR)

#'#Fit a VECM
vecm <- VECM(zeroyld, lag=1, r=1, estim="ML")
logLik(vecm)

Description

Logistic Smooth Transition AutoRegressive model.

Usage

lstar(x, m, d=1, steps=d, series, mL, mH, mTh, thDelay, 
          thVar, th, gamma, trace=TRUE, include = c("const", "trend","none", "both"), 
          control=list(), starting.control=list())

Arguments

Details

$x_{t+s} = ( \phi_{1,0} + \phi_{1,1} x_t + \phi_{1,2} x_{t-d} + \dots + \phi_{1,mL} x_{t - (mL-1)d} ) G( z_t, th, \gamma ) + ( \phi_{2,0} + \phi_{2,1} x_t + \phi_{2,2} x_{t-d} + \dots + \phi_{2,mH} x_{t - (mH-1)d} ) (1 - G( z_t, th, \gamma ) ) + \epsilon_{t+steps}$

with z the threshold variable, and $G$ the logistic function, computed as plogis(q, location = th, scale = 1/gamma), so see plogis documentation for details on the logistic function formulation and parameters meanings. The threshold variable can alternatively be specified by:

mTh

$z[t] = x[t] mTh[1] + x[t-d] mTh[2] + \dots + x[t-(m-1)d] mTh[m]$

thDelay

$z[t] = x[t - thDelay*d ]$

thVar

$z[t] = thVar[t]$

Note that if starting values for phi1 and phi2 are provided, isn't necessary to specify mL and mH. Further, the user has to specify only one parameter between mTh, thDelay and thVar for indicating the threshold variable.

Estimation of the transition parameters th and gamma, as well as the regression parameters phi1 and phi2, is done using concentrated least squares, as suggested in Leybourne et al. (1996).

Given th and gamma, the model is linear, so regression coefficients can be obtained as usual by OLS. So the nonlinear numerical search needs only to be done for th and gamma; the regression parameters are then recovered by OLS again from the optimal th and gamma.

For the nonlinear estimation of the parameters th and gamma, the program uses the optim function, with optimization method BFGS using the analytical gradient. For the estimation of standard values, optim is re-run using the complete Least Squares objective function, and the standard errors are obtained by inverting the hessian. You can pass further arguments to optim directly with the control list argument. For instance, the option maxit maybe useful when there are convergence issues (see examples).

Starting parameters are obtained doing a simple two-dimensional grid-search over th and gamma. Parameters of the grid (interval for the values, dimension of the grid) can be passed to starting.control.

nTh

The number of threshold values (th) in the grid. Defaults to 200

nGamma

The number of smoothing values (gamma) in the grid. Defaults to 40

trim

The minimal percentage of observations in each regime. Defaults to 10% (possible threshold values are between the 0.1 and 0.9 quantile)

gammaInt

The lower and higher smoothing values of the grid. Defaults to c(1,40)

thInt

The lower and higher threshold values of the grid. When not specified (default, i.e NA), the interval are the trim quantiles above.

Value

An object of class nlar, subclass lstar, i.e. a list with fitted model informations.

Author(s)

Antonio, Fabio Di Narzo

References

Non-linear time series models in empirical finance, Philip Hans Franses and Dick van Dijk, Cambridge: Cambridge University Press (2000).

Non-Linear Time Series: A Dynamical Systems Approach, Tong, H., Oxford: Oxford University Press (1990).

Leybourne, S., Newbold, P., Vougas, D. (1998) Unit roots and smooth transitions, Journal of Time Series Analysis, 19: 83-97

See Also

plot.lstar for details on plots produced for this model from the plot generic.

Examples

#fit a LSTAR model. Note 'maxit': slow convergence
mod.lstar <- lstar(log10(lynx), m=2, mTh=c(0,1), control=list(maxit=3000))
mod.lstar

#fit a LSTAR model without a constant in both regimes. 
mod.lstar2 <- lstar(log10(lynx), m=1,  include="none")
mod.lstar2

#Note in example below that the initial grid search seems to be to narrow. 
# Extend it, and evaluate more values (slow!):
controls <- list(gammaInt=c(1,2000), nGamma=50)
mod.lstar3 <- lstar(log10(lynx), m=1,  include="none", starting.control=controls)
mod.lstar3

# a few methods for lstar:
summary(mod.lstar)
residuals(mod.lstar)
AIC(mod.lstar)
BIC(mod.lstar)
plot(mod.lstar)
predict(mod.lstar, n.ahead=5)

Description

Monthly US civilian unemployment rate from 1948 through 2004. Used in the book financial time series for Tsay (2005, chapter 4)

Usage

data(m.unrate)

Format

zoo object

Source

‘⁠https://faculty.chicagobooth.edu/ruey-s-tsay/research/analysis-of-financial-time-series-2nd-edition⁠’

References

Ruey Tsay (2005) Analysis of Financial Time Series, 2nd ed. (Wiley, ch. 4)

Description

This optional function allows the user to set different restrictions for the threshold grid search in function selectSETAR.

Usage

MakeThSpec(
  ngrid = c("All", "Half", "Third", "Quarter"),
  exact = NULL,
  int = c("from", "to"),
  around = "val",
  ...
)

Arguments

Details

This function is just to check the inputs for the specification of the grid search. If not provided, the search will be in the biggest interval (⁠ngrid ="All"⁠) between the minimum and maximum values. The user can reduce it by giving setting "Half" (only every two points is taken) and so on, or setting a number.

The search can also be made around a point, or between two points. When between a point, the argument ngrid is still used, whereas for around, a value of 30 is taken as default value if ngrid is not specified by user.

Value

The input values are given as output after checking for consistency (only one of exact/int/around should be given).

Author(s)

Matthieu Stigler

See Also

selectSETAR

Examples

sun<-(sqrt(sunspot.year+1)-1)*2		
selectSETAR(sun, m=3, th=MakeThSpec(exact=10.40967),criterion="SSR", d=1, thDelay=0:2,
           plot=FALSE, nthresh=1)
#when pre-sepcified value does not correspond, function will search nearest value
selectSETAR(sun, m=3, th=MakeThSpec(exact=10.4),criterion="SSR", d=1, thDelay=0:2,
           plot=FALSE, nthresh=1)
#search around:
selectSETAR(sun, m=3, th=MakeThSpec(around=10.40967, ngrid=20),criterion="SSR", d=1, thDelay=0:2,
           plot=FALSE, nthresh=1)
#search in an interval
selectSETAR(sun, m=3, th=MakeThSpec(int=c(10, 11), ngrid=20),criterion="SSR", d=1, thDelay=0:2,
           plot=FALSE, nthresh=1)
#reduce size of the grid:
selectSETAR(sun, m=3, th=MakeThSpec(ngrid="Half"),criterion="SSR", d=1, thDelay=0:2,
           plot=FALSE, nthresh=1)


# 2 thresholds:
selectSETAR(sun, m=3, th=MakeThSpec(ngrid="Half"),criterion="SSR", d=1, thDelay=0:2,
           plot=FALSE, nthresh=2)

Description

Generic function to compute the Mean Absolute Percent Error of a fitted model.

Usage

MAPE(object, ...)

## Default S3 method:
MAPE(object, ...)

Arguments

Value

Computed Mean Absolute Percent Error for the fitted model.

Author(s)

Antonio, Fabio Di Narzo

Description

Generic function to compute the Mean Squared Error of a fitted model.

Usage

mse(object, ...)

## Default S3 method:
mse(object, ...)

Arguments

Value

Computed MSE for the fitted model.

Author(s)

Antonio, Fabio Di Narzo

Description

Generic ‘nlar’ methods. Method ‘nlar’ is described in a separate page: nlar

Usage

## S3 method for class 'nlar'
coef(object, ...)

## S3 method for class 'nlar'
fitted(object, ...)

## S3 method for class 'nlar'
residuals(object, initVal = TRUE, timeAttr = TRUE, ...)

## S3 method for class 'nlar'
deviance(object, ...)

## S3 method for class 'nlar'
mse(object, ...)

## S3 method for class 'nlar'
AIC(object, k = 2, ...)

## S3 method for class 'nlar'
BIC(object, ...)

## S3 method for class 'nlar'
MAPE(object, ...)

## S3 method for class 'nlar'
summary(object, ...)

## S3 method for class 'nlar'
plot(x, ask = interactive(), ...)

## S3 method for class 'nlar'
toLatex(object, digits, label, ...)

Arguments

Details

MAPE

Mean Absolute Percent Error

mse

Mean Square Error

plot

Diagnostic plots

Author(s)

Antonio, Fabio Di Narzo

See Also

availableModels for listing all currently available models.

Examples

x <- log10(lynx)
mod.setar <- setar(x, m=2, thDelay=1, th=3.25)
mod.setar
AIC(mod.setar)
mse(mod.setar)
MAPE(mod.setar)
coef(mod.setar)
summary(mod.setar)

e <- residuals(mod.setar)
e <- e[!is.na(e)]
plot(e)
acf(e)

plot(x)
lines(fitted(mod.setar), lty=2)
legend(x=1910, y=3.9,lty=c(1,2), legend=c("observed","fitted"))

plot(mod.setar)

Description

Neural Network nonlinear autoregressive model.

Usage

nnetTs(x, m, d = 1, steps = d, series, size, 
	control = list(trace = FALSE))

Arguments

Details

Neural network model with 1 hidden layer and linear output:

$x_{t+s} = \beta_0 + \sum_{j=1}^D \beta_j g( \gamma_{0j} + \sum_{i=1}^{m} \gamma_{ij} x_{t-(i-1) d} )$

Model is estimated using the nnet function in nnet package. Optimization is done via the BFGS method of optim. Note that for this model, no additional model-specific summary and plot methods are made available from this package.

Value

An object of class nlar, subclass nnetTs, i.e. a list with mostly nnet::nnet internal structures.

Author(s)

Antonio, Fabio Di Narzo

References

Non-linear time series models in empirical finance, Philip Hans Franses and Dick van Dijk, Cambridge: Cambridge University Press (2000).

Non-Linear Time Series: A Dynamical Systems Approach, Tong, H., Oxford: Oxford University Press (1990).

Chaos: A Statistical Perspective, Chan, K., Tong, H., New York: Springer-Verlag (2001).

Examples

#fit a Neural Network model
mod.nnet <- nnetTs(log(lynx), m=2, size=3)
mod.nnet

Description

Plotting methods ‘setar’ and ‘lstar’ subclasses

Usage

## S3 method for class 'setar'
plot(x, ask=interactive(), legend=FALSE, regSwStart, regSwStop, ...)
## S3 method for class 'lstar'
plot(x, ask=interactive(), legend=FALSE, regSwStart, regSwStop, ...)

Arguments

Details

These plot methods produce a plot which gives to you an idea of the behaviour of the fitted model.

Firstly, if embedding dimension is, say, m, m scatterplots are produced. On the x axis you have the lagged time series values. On the y axis the ‘response’ time series values. Observed points are represented with different colors-symbols depending on the level of the threshold variable. Specifically, for the setar model, black means ‘low regime’, red means ‘high regime’. For the lstar model, where the self-threshold variable is continuous, threshold values are grouped in 5 different zones with the same number of points in each. Note that if more than 300 points are to be plotted, they all share the same symbol, and regimes can be distinguished only by color. If you want, by specifying legend=TRUE a legend is added at the upper-left corner of each scatterplot. To each scatterplot, a dashed line is superposed, which links subsequent fitted values.

Finally, a new time series plot is produced, with lines segments coloured depending on the regime (colors meanings are the same of those in the preceding scatterplots). Optionally, you can specify a starting and ending time indices, for zooming on a particular segment of the time series.

Author(s)

Antonio, Fabio Di Narzo

See Also

setar, lstar

nlar-methods for other generic available methods for this kind of objects.

Examples

##
##See 'setar' examples
##

Description

This plot shows how variables in a (T)VECM respond to deviations from the long-term equilibrium

Usage

plot_ECT(x, add.legend = TRUE, legend.location = "topright", ...)

Arguments

Value

a plot, and invisibly the underlying data.frame, containing the ECT and the response for each variable

Examples

data(zeroyld)
vec_l1 <- VECM(zeroyldMeta[, c("long.run", "short.run")], lag =1)
tvec_l1 <- TVECM(zeroyldMeta[, c("long.run", "short.run")], lag =1, 
                plot  = FALSE, trace = FALSE, th1 = list(exact = -1.263))

plot_ECT(vec_l1)
plot_ECT(tvec_l1, legend.location = "bottomright")

Description

Forecasts a VAR or VECM by discarding a part of the sample, and recursively generating a series of updated forecasts.

Usage

predict_rolling(object, ...)

## S3 method for class 'nlVar'
predict_rolling(object, nroll = 10, n.ahead = 1, refit.every, newdata, ...)

Arguments

Details

This function allows to check the out-of sample forecasting accuracy by estimating the model on a sub-sample of the original, then making nroll forecasts of horizon n.ahead, each time by updating the sample. In other words, with a given model estimated on 100 observations, the function will estimate it on say 90 first obs (nroll=10), generate a say 1 step-ahead (n.ahead=1) from obs 90, then using true value 91 for to predict value at 92, etc, till the full sample is used.

Unlike with usual predict() methods, specifiying n.ahead=2 will not generate a 1 step-ahead and a 2 step-ahead forecasts, but only nroll 2 step-ahead forecasts.

Note that while the forecasts are updated with new values, the model estimation is (by default) not updated. This can however be done with the argument fit.every, specifiying at which frequency the model should be re-estimated. By setting it to 1 for example, each time a new observation is taken, the model is reestimated. This is similar to the ugarchroll in package rugarch.

Value

A matrix containing the forecasts.

Author(s)

Matthieu Stigler

See Also

predict.nlar for the standard predict function.

Examples

data(barry)

## model estimated on full sample:
mod_vec <- VECM(barry, lag=2)

## generate 10 1-step-ahead forecasts:
preds_roll <- predict_rolling(mod_vec, nroll=10)

## plot the results:
plot(window(barry[,"dolcan"],start=1998), type="l", ylab="barry: dolcan")
preds_roll_ts <- ts(preds_roll$pred, start=time(barry)[nrow(barry)-10], freq=12)
lines(preds_roll_ts[,"dolcan"], col=2, lty=2)
legend("bottomright", lty=c(1,2), col=1:2, leg=c("True", "Fitted"))
title("Comparison of true and rolling 1-ahead forecasts\n")

Description

Forecasting a non-linear model object of general class ‘nlar’, including ‘setar’ and ‘star’.

Usage

## S3 method for class 'nlar'
predict(
  object,
  newdata,
  n.ahead = 1,
  type = c("naive", "MC", "bootstrap", "block-bootstrap"),
  nboot = 100,
  ci = 0.95,
  block.size = 3,
  boot1Zero = TRUE,
  seed = NULL,
  ...
)

Arguments

Details

The forecasts are obtained recursively from the estimated model. Given that the models are non-linear, ignoring the residuals in the 2- and more steps ahead forecasts leads to biased forecasts (so-called naive). Different resampling methods, averaging n.boot times over future residuals, are available:

naive

No residuals

MC

Monte-Carlo method, where residuals are taken from a normal distribution, with a standard deviation equal to the residuals sd.

bootstrap

Residuals are resampled from the empirical residuals from the model.

block-bootstrap

Same as bootstrap, but residuals are resampled in block, with size block.size

The MC and bootstrap methods correspond to equations 3.90 and 3.91 of Franses and van Dijk (2000, p. 121). The bootstrap/MC is initiated either from the first forecast, n.ahead=1 (set with boot1zero to TRUE), or from the second only.

When the forecast method is based on resampling, forecast intervals are available. These are obtained simply as empirical ci quantiles of the resampled forecasts (cf Method 2 in Franses and van Dijk, 2000, p. 122).

Value

A ‘ts’ object, or, in the case of MC/bootstrap, a list containing the prediction (pred) and the forecast standard errors (se).

Author(s)

Matthieu Stigler

References

Non-linear time series models in empirical finance, Philip Hans Franses and Dick van Dijk, Cambridge: Cambridge University Press (2000).

See Also

The model fitting functions setar, lstar.

A more sophisticated predict function, allowing to do sub-sample rolling predictions: predict_rolling.

Examples

x.train <- window(log10(lynx), end = 1924)
x.test <- window(log10(lynx), start = 1925)


### Use different forecasting methods:
mod.set <- setar(x.train, m=2, thDelay=0)
pred_setar_naive <- predict(mod.set, n.ahead=10)
pred_setar_boot <- predict(mod.set, n.ahead=10, type="bootstrap", n.boot=200)
pred_setar_Bboot <- predict(mod.set, n.ahead=10, type="block-bootstrap", n.boot=200)
pred_setar_MC <- predict(mod.set, n.ahead=10, type="bootstrap", n.boot=200)

## Plot to compare results:
pred_range <- range(pred_setar_naive, pred_setar_boot$pred, pred_setar_MC$pred, na.rm=TRUE)
plot(x.test, ylim=pred_range, main="Comparison of forecasts methods from same SETAR")
lines(pred_setar_naive, lty=2, col=2)
lines(pred_setar_boot$pred, lty=3, col=3)
lines(pred_setar_Bboot$pred, lty=4, col=4)
lines(pred_setar_MC$pred, lty=5, col=5)

legLabels <- c("Observed", "Naive F", "Bootstrap F","Block-Bootstrap F", "MC F")
legend("bottomleft", leg=legLabels, lty=1:5, col=1:5)

Description

Forecasting the level of a series estimated by ‘VAR’ / ‘VECM’ or ‘TVAR’

Usage

## S3 method for class 'TVAR'
predict(object, newdata, n.ahead = 5, newdataTrendStart, ...)

## S3 method for class 'VAR'
predict(object, newdata, n.ahead = 5, newdataTrendStart, exoPred = NULL, ...)

Arguments

Details

The forecasts are obtained recursively, and are for the levels of the series.

When providing newdata, newdata has to be ordered chronologically, so that the first row/element is the earliest value.

For VECM, the forecasts are obtained by transforming the VECM to a VAR (using function VARrep). Note that a VECM(lag=p) corresponds to a VAR(lag=p+1), so that if the user provides newdata for a VECM(lag=p), newdata should actually contain p+1 rows.

Value

A matrix of predicted values.

Author(s)

Matthieu Stigler

See Also

lineVar and VECM. VARrep

A more sophisticated predict function, allowing to do sub-sample rolling predictions: predict_rolling.

Examples

data(barry)
barry_in <- head(barry, -5)
barry_out <- tail(barry, 5)

mod_vecm <- VECM(barry_in, lag=2)
mod_var <- lineVar(barry_in, lag=3)
mod_tvar <- TVAR(barry_in, lag=3, nthresh=1, thDelay=1)

pred_vecm <- predict(mod_vecm)
pred_var  <- predict(mod_var) 
pred_tvar <- predict(mod_tvar)

 
## compare forecasts on a plot
n <- 30
plot(1:n, tail(barry[,1], n), type="l", xlim=c(0,n))
lines((n-5+1):n, pred_var[,1], lty=2, col=2)
lines((n-5+1):n, pred_vecm[,1], lty=2, col=3)
lines((n-5+1):n, pred_tvar[,1], lty=2, col=4) 
legend("bottomright", lty=c(1,2,2,2), col=1:4, legend=c("true", "var", "vecm", "tvar"))

## example for newdata:
all.equal(predict(mod_vecm), predict(mod_vecm, newdata=barry[c(317, 318, 319),]))

Description

Selection of the cointegrating rank and the lags with Information criterion (AIC, BIC).

Usage

rank.select(
  data,
  lag.max = 10,
  r.max = ncol(data) - 1,
  include = c("const", "trend", "none", "both"),
  fitMeasure = c("SSR", "LL"),
  sameSample = TRUE,
  returnModels = FALSE
)

## S3 method for class 'rank.select'
print(x, ...)

## S3 method for class 'rank.select'
as.data.frame(x, ...)

## S3 method for class 'rank.select'
summary(object, ...)

Arguments

Details

This function estimates the AIC, BIC and Hannan-Quinn for each rank (up to lags.max) and lags (up to lags.max). This method has been shown to be useful to select simultaneously the rank and the lags, see references.

Value

An object of class ‘rank.select’, with ‘print’ and ‘summary methods’, containing among other the matrices of AIC/BIC/HQ, the Likelihood, and best ranks according to each criterion.

Author(s)

Matthieu Stigler

References

- Aznar A and Salvador M (2002). Selecting The Rank Of The Cointegration Space And The Form Of The Intercept Using An Information Criterion. Econometric Theory, *18*(04), pp. 926-947. <URL: http://ideas.repec.org/a/cup/etheor/v18y2002i04p926-947_18.html>.

-Cheng X and Phillips PCB (2009). Semiparametric cointegrating rank selection. Econometrics Journal , *12*(s1), pp. S83-S104. <URL: http://ideas.repec.org/a/ect/emjrnl/v12y2009is1ps83-s104.html>.

- Gonzalo J and Pitarakis J (1998). Specification via model selection in vector error correction models. Economics Letters, *60*(3), pp. 321 - 328. ISSN 0165-1765, <URL: http://dx.doi.org/DOI: 10.1016/S0165-1765(98)00129-3>.

- Kapetanios G (2004). The Asymptotic Distribution Of The Cointegration Rank Estimator Under The Akaike Information Criterion. Econometric Theory, *20*(04), pp. 735-742. <URL: http://ideas.repec.org/a/cup/etheor/v20y2004i04p735-742_20.html>.

- Wang Z and Bessler DA (2005). A Monte Carlo Study On The Selection Of Cointegrating Rank Using Information Criteria. Econometric Theory, *21*(03), pp. 593-620. <URL: http://ideas.repec.org/a/cup/etheor/v21y2005i03p593-620_05.html>.

See Also

VECM for estimating a VECM. rank.test (or ca.jo in package urca) for the classical Johansen cointegration test.

Examples

data(barry)

# 
rk_sel <- rank.select(barry)
rk_sel
summary(rk_sel)

Description

Maximum-likelihood test of the cointegrating rank.

Usage

rank.test(vecm, type = c("eigen", "trace"), r_null, cval = 0.05)

## S3 method for class 'rank.test'
print(x, ...)

## S3 method for class 'rank.test'
summary(object, digits = max(1, getOption("digits") - 3), ...)

Arguments

Details

This function computes the two maximum-likelihood tests for the cointegration rank from Johansen (1996). Tests are:

trace

Test the hypothesis of rank ‘h’ against rank ‘K’, i.e. against the alternative that the system is stationary.

eigenvalue

Test the hypothesis of rank ‘h’ against rank ‘h+1’.

The test works for five specifications of the deterministic terms as in Doornik et al (1998), to be specified in the previous call to VECM:

H_ql

Unrestricted constant and trend: use include="both"

H_l

Unrestricted constant and restricted trend: use include="const"

and LRinclude="trend"

H_lc

Unrestricted constant and no trend: use include="const"

H_c

Restricted constant and no trend: use LRinclude="const"

H_z

No constant nor trend: use include="none"

Two testing procedures can be used:

Specific test

By specifying a value for ‘r_null’. The ‘pval’ value returned gives the speciifc p-value.

Automatic test

If not value is specified for ‘r_null’, the function makes a simple automatic test: returns the rank (slot ‘r’) of the first test not rejected (level specified by arg cval) as recommend i.a. in Doornik et al (1998, p. 544).

A full table with both test statistics ad their respective p-values is given in the summary method.

P-values are obtained from the gamma approximation from Doornik (1998, 1999). Small sample values adjusted for the sample site are also available in the summary method. Note that the ‘effective sample size’ for the these values is different from output in gretl for example.

Value

An object of class ‘rank.test’, with ‘print’ and ‘summary methods’.

Comparison with urca

While ca.jo in package urca and rank.test both implement Johansen tests, there are a few differences:

• rank.test gives p-values, while ca.jo gives only critical values.

• rank.test allows for five different specifications of deterministic terms (see above), ca.jo for only three.

• ca.jo allows for seasonal and exogenous regressors, which is not available in rank.test.

• The lag is specified differently: K from ca.jo corresponds to lag+1 in rank.test.

Author(s)

Matthieu Stigler

References

- Doornik, J. A. (1998) Approximations to the Asymptotic Distributions of Cointegration Tests, Journal of Economic Surveys, 12, 573-93

- Doornik, J. A. (1999) Erratum [Approximations to the Asymptotic Distribution of Cointegration Tests], Journal of Economic Surveys, 13, i

- Doornik, Hendry and Nielsen (1998) Inference in Cointegrating Models: UK M1 Revisited, Journal of Economic Surveys, 12, 533-72

- Johansen, S. (1996) Likelihood-based inference in cointegrated Vector Autoregressive Models, Oxford University Press

See Also

VECM for estimating a VECM. rank.select to estimate the rank based on information criteria.

ca.jo in package urca for another implementation of Johansen cointegration test (see section ‘Comparison with urca’ for more infos).

Examples

data(barry)

## estimate the VECM with Johansen! 
ve <- VECM(barry, lag=1, estim="ML")

## specific test:
ve_test_spec <- rank.test(ve, r_null=1)
ve_test_spec_tr <- rank.test(ve, r_null=1, type="trace")

ve_test_spec
ve_test_spec_tr

## No specific test: automatic method
ve_test_unspec <- rank.test(ve)
ve_test_unspec_tr <- rank.test(ve, type="trace")

ve_test_unspec
ve_test_unspec_tr

## summary method: output will be same for all types/test procedure:
summary(ve_test_unspec_tr)

## The function works for many specification of the VECM(), try:
rank.test(VECM(barry, lag=3, estim="ML"))
rank.test(VECM(barry, lag=3, include="both",estim="ML"))
rank.test(VECM(barry, lag=3, LRinclude="const",estim="ML"))

## Note that the tests are simple likelihood ratio, and hence can be obtained also manually:
-2*(logLik(ve, r=1)-logLik(ve, r=2)) # eigen test, 1 against 2
-2*(logLik(ve, r=1)-logLik(ve, r=3)) # eigen test, 1 against 3

Description

This function allows to extract the indicator variable specifying the regime in which the process is at time t.

Usage

regime(object, initVal = TRUE, timeAttr = TRUE, series = NULL, ...)

## S3 method for class 'lstar'
regime(object, initVal = TRUE, timeAttr = TRUE, series, discretize = TRUE, ...)

Arguments

Value

Time series of same attributes as input to setar.

Author(s)

Matthieu Stigler

Examples

set<-setar(lynx, m=3)
regime(set)
regime(set, time=FALSE, initVal=FALSE)

plot(regime(set))

Description

Bootstrap a vector according to multiple resampling schemes: resampling, block resampling, Wild bootstrap.

Usage

resample_vec(
  x,
  boot.scheme = c("resample", "resample_block", "wild1", "wild2", "check"),
  seed = NULL,
  block.size = 2
)

Arguments

Details

This function offers various bootstrap/resampling schemes:

resample

Resampling with replacement

resample_block

Resampling contiguous observations (blocks) with replacement. Use argument block.size

wild1

Wild bootstrap: do not resample, but add a N(0,1) distribution to each value

wild12

Wild bootstrap: same, but add instead -1 or 1.

Description

Extracts the global and regime-dependent variance of the residuals

Usage

resVar(x, adj=c("OLS", "ML"))

Arguments

Details

The degree of freedom adjustment in the formula for the variance is the number of parameters when adj="OLS" or zero when adj="ML".

Value

A vector containing:

Author(s)

Matthieu Stigler

References

Non-Linear Time Series: A Dynamical Systems Approach, Tong, H., Oxford: Oxford University Press (1990).

Examples

#Lynx model as in Tong (1980, p. 387)
mod.setar <- setar(log10(lynx), mL=7,mH=2, thDelay=1, th=3.116)
summary(mod.setar)
#coefficients are same for lower regime but differ for higer

resVar(mod.setar, adj="ML")
#variance or the residuals is same for lower regime but differ for higer regime and hence for total

#Lynx model as in Tong (1980, p. 405)
mod.setar2 <- setar(log10(lynx), mL=1,mM=7,mH=2, thDelay=1, nthresh=2,th=c(2.373, 3.154))
round(coefficients(mod.setar2),3)

resVar(mod.setar2, adj="ML")

Description

Automatic selection of model hyper-parameters

Usage

selectLSTAR(x, m, d=1, steps=d, mL = 1:m, mH = 1:m, thDelay=0:(m-1), 
            fast=TRUE, trace=FALSE)
selectNNET(x, m, d=1, steps=d, size=1:(m+1), maxit=1e3, trace=FALSE)

Arguments

Details

Functions for automatic selection of LSTAR and NNET models hyper parameters. An exhaustive search over all possible combinations of values of specified hyper-parameters is performed. Embedding parameters m,d,steps are kept fixed.

Selection criterion is the usual AIC.

For the LSTAR model, two methods are offered:

fast=FALSE

Each model is run separately, each time using the full grid search for starting values.

fast=TRUE

Only the first model is run with a full grid search, while the subsequent use the first model results for their starting values.

Value

A data-frame, with columns giving hyper-parameter values and the computed AIC for each row (only the best 10s are returned)

Author(s)

Antonio, Fabio Di Narzo

Examples

llynx <- log10(lynx)
selectLSTAR(llynx, m=2)
selectNNET(llynx, m=3, size=1:5)

Description

Automatic selection of SETAR hyper-parameters

Usage

selectSETAR(x, m, d=1, steps=d, series, mL, mH,mM, thDelay=0, mTh, thVar, 
	    th=MakeThSpec(), trace=TRUE, include = c("const", "trend","none", "both"), 
	    common=c("none", "include","lags", "both"), model=c("TAR", "MTAR"), 
	    ML=seq_len(mL), MH=seq_len(mH), MM=seq_len(mM), nthresh=1, trim=0.15,
	    criterion = c("pooled-AIC", "AIC","BIC", "SSR"),  plot=TRUE,
	    max.iter=2, type=c("level", "diff", "ADF"), same.lags=FALSE, 
	    restriction=c("none","OuterSymAll","OuterSymTh"),  hpc=c("none", "foreach") )

Arguments