Nothing
#' VLSTAR- Estimation
#'
#' This function allows the user to estimate the coefficients of a VLSTAR model with \emph{m} regimes through maximum likelihood or nonlinear least squares.
#' The set of starting values of Gamma and C for the convergence algorithm can be either passed or obtained via searching grid.
#'
#' The multivariate smooth transition model is an extension of the smooth transition regression model introduced by Bacon and Watts (1971) (see also Anderson and Vahid, 1998). The general model is
#' \deqn{y_{t} = \mu_0+\sum_{j=1}^{p}\Phi_{0,j}\,y_{t-j}+A_0 x_t \cdot G_t(s_t;\gamma,c)[\mu_{1}+\sum_{j=1}^{p}\Phi_{1,j}\,y_{t-j}+A_1x_t]+\varepsilon_t}
#' where \eqn{\mu_{0}} and \eqn{\mu_{1}} are the \eqn{\tilde{n} \times 1} vectors of intercepts, \eqn{\Phi_{0,j}} and \eqn{\Phi_{1,j}} are square
#' \eqn{\tilde{n}\times\tilde{n}} matrices of parameters for lags \eqn{j=1,2,\dots,p}, A_0 and A_1 are \eqn{\tilde{n}\times k} matrices of parameters,
#' x_t is the \eqn{k \times 1} vector of exogenous variables and \eqn{\varepsilon_t} is the innovation. Finally, \eqn{G_t(s_t;\gamma,c)} is a \eqn{\tilde{n}\times \tilde{n}} diagonal matrix of transition function at time \emph{t}, such that
#' \deqn{G_t(s_t;\gamma,c)=\{G_{1,t}(s_{1,t};\gamma_{1},c_{1}),G_{2,t}(s_{2,t};\gamma_{2},c_{2}),
#' \dots,G_{\tilde{n},t}(s_{\tilde{n},t};\gamma_{\tilde{n}},c_{\tilde{n}})\}.}
#'
#' Each diagonal element \eqn{G_{i,t}^r} is specified as a logistic cumulative density functions, i.e.
#' \deqn{G_{i,t}^r(s_{i,t}^r; \gamma_i^r, c_i^r) = \left[1 + \exp\big\{-\gamma_i^r(s_{i,t}^r-c_i^r)\big\}\right]^{-1}}
#' for \eqn{latex}{i = 1,2, \dots, \tilde{n}} and \eqn{r=0,1,\dots,m-1}, so that the first model is a Vector Logistic Smooth Transition AutoRegressive
#' (VLSTAR) model.
#' The ML estimator of \eqn{\theta} is obtained by solving the optimization problem
#' \deqn{\hat{\theta}_{ML} = arg \max_{\theta}log L(\theta)}
#' where \eqn{log L(\theta)} is the log-likelihood function of VLSTAR model, given by
#' \deqn{ ll(y_t|I_t;\theta)=-\frac{T\tilde{n}}{2}\ln(2\pi)-\frac{T}{2}\ln|\Omega|-\frac{1}{2}\sum_{t=1}^{T}(y_t-\tilde{G}_tB\,z_t)'\Omega^{-1}(y_t-\tilde{G}_tB\,z_t)}
#' The NLS estimators of the VLSTAR model are obtained by solving the optimization problem
#' \deqn{\hat{\theta}_{NLS} = arg \min_{\theta}\sum_{t=1}^{T}(y_t - \Psi_t'B'x_t)'(y_t - \Psi_t'B'x_t).}
#' Generally, the optimization algorithm may converge to some local minimum. For this reason, providing valid starting values of \eqn{\theta} is crucial. If there is no clear indication on the initial set of parameters, \eqn{\theta}, this can be done by implementing a grid search. Thus, a discrete grid in the parameter space of \eqn{\Gamma} and C is create to obtain the estimates of B conditionally on each point in the grid. The initial pair of \eqn{\Gamma} and C producing the smallest sum of squared residuals is chosen as initial values, then the model is linear in parameters.
#' The algorithm is the following:
#' \enumerate{
#' \item Construction of the grid for \eqn{\Gamma} and C, computing \eqn{\Psi} for each poin in the grid
#' \item Estimation of \eqn{\hat{B}} in each equation, calculating the residual sum of squares, \eqn{Q_t}
#' \item Finding the pair of \eqn{\Gamma} and C providing the smallest \eqn{Q_t}
#' \item Once obtained the starting-values, estimation of parameters, \emph{B}, via nonlinear least squares (NLS)
#' \item Estimation of \eqn{\Gamma} and C given the parameters found in step 4
#' \item Repeat step 4 and 5 until convergence.
#' }
#'
#' @param y \code{data.frame} or \code{matrix} of dependent variables of dimension \code{(Txn)}
#' @param exo (optional) \code{data.frame} or \code{matrix} of exogenous variables of dimension \code{(Txk)}
#' @param p lag order
#' @param m number of regimes
#' @param st single transition variable for all the equation of dimension \code{(Tx1)}
#' @param constant \code{TRUE} or \code{FALSE} to include or not the constant
#' @param starting set of intial values for Gamma and C, inserted as a list of length \code{m-1}.
#' Each element of the list should contain a \code{data.frame} with 2 columns (one for Gamma and one for c), and \code{n} rows.
#' @param method Fitting method: maximum likelihood or nonlinear least squares.
#' @param singlecgamma \code{TRUE} or \code{FALSE} to use single gamma and c
#' @param n.iter number of iteration of the algorithm until forced convergence
#' @param epsilon convergence check measure
#' @param ncores Number of cores used for parallel computation. Set to \code{NULL} by default and automatically calculated.
#' @return An object of class \code{VLSTAR}, with standard methods.
#' @references Anderson H.M. and Vahid F. (1998), Testing multiple equation systems for common nonlinear components. \emph{Journal of Econometrics}. 84: 1-36
#'
#' Bacon D.W. and Watts D.G. (1971), Estimating the transition between two intersecting straight lines. \emph{Biometrika}. 58: 525-534
#'
#' Terasvirta T. and Yang Y. (2014), Specification, Estimation and Evaluation of Vector Smooth Transition Autoregressive Models with Applications. \emph{CREATES Research Paper 2014-8}
#' @author Andrea Bucci
#' @export
#' @keywords VLSTAR
#' @examples
#' \donttest{
#' data(Realized)
#' y <- Realized[-1,1:10]
#' y <- y[-nrow(y),]
#' st <- Realized[-nrow(Realized),1]
#' st <- st[-length(st)]
#' stvalues <- startingVLSTAR(y, p = 1, n.combi = 3,
#' singlecgamma = FALSE, st = st, ncores = 1)
#'fit.VLSTAR <- VLSTAR(y, p = 1, singlecgamma = FALSE, starting = stvalues,
#' n.iter = 1, st = st, method ='NLS', ncores = 1)
#'# a few methods for VLSTAR
#'print(fit.VLSTAR)
#'summary(fit.VLSTAR)
#'plot(fit.VLSTAR)
#'predict(fit.VLSTAR, st.num = 1, n.ahead = 1)
#'logLik(fit.VLSTAR, type = 'Univariate')
#'coef(fit.VLSTAR)}
#'
VLSTAR <- function(y, exo = NULL, p = 1,
m = 2, st = NULL, constant = TRUE, starting = NULL,
method = c('ML', 'NLS'), n.iter = 500,
singlecgamma = FALSE,
epsilon = 10^(-3), ncores = NULL){
chk <- Sys.getenv("_R_CHECK_LIMIT_CORES_", "")
if(is.null(ncores)){
if (nzchar(chk) && chk == "TRUE") {
# use 2 cores in CRAN/Travis/AppVeyor
ncores <- 2L
} else {
ncores <- parallel::detectCores()
}}
y <- as.matrix(y)
x <- exo
method <- match.arg(method)
##Checks and warnings
if (anyNA(y))
stop("\nNAs in y.\n")
if(m < 2)
stop('The number of regimes should be greater than one.')
if(is.null(st))
stop('The transition variable must be supplied.')
if(is.null(x)){
if(length(y[,1]) != length(st))
stop('The length of the variables does not match!')
}else{
if(length(y[,1]) != length(as.matrix(x[,1])) | length(st) != length(as.matrix(x[,1])) | length(y[,1]) != length(st))
stop('The length of the variables does not match!')
}
if(is.null(starting)){
stop('Starting values should be provided.')
}
if(!is.list(starting)){
stop('Starting c and gamma should be put in a list.')
}
if(is.null(p) || p < 1){
stop('Please, specify a valid lag order.')
}
if (ncol(y) < 2)
stop("The matrix 'y' should contain at least two variables. For univariate analysis consider lstar() function in this package.\n")
if (is.null(colnames(y))) {
colnames(y) <- paste("y", 1:ncol(y), sep = "")
warning(paste("No column names supplied in y, using:",
paste(colnames(y), collapse = ", "), ", instead.\n"))
}
colnames(y) <- make.names(colnames(y))
##Definition of dimensions, creating variable x with constant
yt <- zoo::zoo(y)
ylag <- stats::lag(yt, -(1:p))
ylag <- as.matrix(ylag)
if(p>1){
lagg <- p-1
ylag <- ylag[-(1:lagg),]
}
y <- y[-c(1:p), ]
ncoly <- ncol(y)
if(singlecgamma == TRUE){
starting1 <- list()
for(i in 1:(m-1)){
starting1[[i]] <- cbind(rep(starting[[i]][,1], ncoly), rep(starting[[i]][,2], ncoly))
}
starting = starting1
}
if(!is.null(starting)){
if(length(starting)!= (m-1)){
stop('The length of the list of initial values should be equal to m-1.')
}else{
if(any(unlist(lapply(starting, ncol))!=2) | any(unlist(lapply(starting, nrow))!=ncoly)){
stop('Each element of the starting argument should have two columns and n rows.')
}
}
}
ncolylag <- ncoly*p
nrowy <- nrow(y)
ncolx1 <- ncol(x)
const <- rep(1, nrowy)
if (constant == TRUE){
if(!is.null(exo)){
x1a <- as.matrix(x[-c(1:p),])
x <- as.matrix(cbind(const,ylag,x1a))
}else{
x <- as.matrix(cbind(const,ylag))
}
ncolx <- ncol(x)
} else{
x <- as.matrix(x)
}
nrowx <- nrow(x)
st <- st[(1+p):length(st)]
ncolx <- ncol(x)
ny <- ifelse(singlecgamma == TRUE, 1, ncoly)
q <- ncol(x)-ncolylag
PARAM <- starting
####Estimating VLSTAR model####
##Estimating initial values to be used in the iterative algorithm
In <- diag(ncoly)
glog <- matrix(ncol=ny, nrow = nrowy)
GT <- list()
Gtilde <- list()
kro <- list()
for (i in 1:nrowx){
for (t in 1:(m-1)){
for (j in 1 : ny){
glog[i,j] <- (1L+exp(-PARAM[[t]][j,1]*(st[i]-PARAM[[t]][j,2])))^(-1)
}
if(singlecgamma == TRUE){
Gt <- diag(rep(glog[i,1], ncoly))
}else{
Gt <- diag(glog[i,])
}
GT[[t]] <- Gt
}
Gtilde[[i]] <- t(cbind(In, do.call(cbind,GT)))
kro[[i]] <- kronecker(Gtilde[[i]], x[i,])
}
M <- t(do.call("cbind", kro))
Y <- matrixcalc::vec(t(y))
Bhat <- MASS::ginv(t(M)%*%M)%*%t(M)%*%Y
BB <- ks::invvec(Bhat, ncol = (m*ncoly), nrow = (ncolylag + q)) ##Estimated coefficients
##Calculating the estimated covariance matrix
resi <- list()
Ehat1 <- matrix(NA, ncol = ncoly, nrow = nrowy)
for (i in 1:nrowx){
resi[[i]] <- y[i, ] - t(Gtilde[[i]])%*%t(BB)%*%x[i,]
}
Ehat1 <- t(do.call("cbind", resi))
#Estimated covariance matrix
Omegahat <- (t(Ehat1)%*%Ehat1)/(nrowy)
data = list(y = y, x = x, m = m, BB = BB, Omegahat = Omegahat, st = st, singlecgamma = singlecgamma)
#Inizialization of iter
iter <- 0
ll0 <- 10^(4)
epsi <- epsilon ##Value used as convergence check
loglik1 <- NULL
omega <- list()
PARAM1 <- list()
for(t in 1:(m-1)){
PARAM1[[t]] <- as.data.frame(PARAM[[t]])
}
param <- do.call(rbind,PARAM1)
param <- matrixcalc::vec(as.matrix(param)) ##the parameters c and gamma are vectorized in order to be passed in the optimParallel function
err <- 10^5
##Actual iteration to estimate the coefficients
# cl <- parallel::makeCluster(ncores, outfile="") # set the number of processor cores
# parallel::setDefaultCluster(cl=cl)
if(method == 'ML'){
#NLS Estimation of Bhat and Omegahat to be used in the first iteration of maximum likelihood
message('Maximum likelihood estimation\n')
errdif <- 10^5
#Convergence algorithm
#1. Maximum likelihood estimation of gamma and c with NLS estimates of Bhat and Omegahat
#2. Maximum likelihood estimation of Bhat with new values of gamma and c
#3. Convergence check
#4. 1-2-3 until convergence
while (iter < n.iter & errdif > epsi){
Sys.sleep(0)
iter <- iter+1
#Parameters
low1 <- replicate(ncoly, 0)
#1.Maximum likelihood estimation of gamma and c
cl <- parallel::makeCluster(ncores) # set the number of processor cores
parallel::setDefaultCluster(cl=cl)
param1 <- optimParallel::optimParallel(par = as.vector(param), fn = loglike, lower = c(low1, apply(y, 2, min)),
data = data, parallel = list(cl = cl, forward = FALSE, loginfo = FALSE))
parallel::stopCluster(cl)
cgam1 <- matrix(param1$par, ncol = 2, nrow = (ny*(m-1)))
#2.Maximum likelihood estimation of Bhat with new values of gamma and c
glog <- rep(0, ncol(y))
GT <- list()
Gtilde <- list()
XX <- list()
XYOG <- list()
kro <- list()
PsiOmegaPsi <- list()
for (i in 1:nrowx){
for(t in 1:(m-1)){
for (j in 1 : ny){
glog[j] <- (1+exp(-cgam1[j,1]*(st[i]-cgam1[j,2])))^(-1)
}
if(singlecgamma == TRUE){
GT[[t]] <- diag(rep(glog[1], ncoly))
}else{
GT[[t]] <- diag(glog)
}
}
Gtilde[[i]] <- t(cbind(In, do.call(cbind,GT)))
XX[[i]] <- x[i,] %*%t(x[i,])
XYOG[[i]] <- matrixcalc::vec((x[i, ]%*%t(y[i,]))%*%MASS::ginv(Omegahat)%*%t(Gtilde[[i]]))
PsiOmegaPsi[[i]] <- Gtilde[[i]]%*%MASS::ginv(Omegahat)%*%t(Gtilde[[i]])
kro[[i]] <- kronecker(PsiOmegaPsi[[i]], XX[[i]])
}
xyog <- Reduce(`+`, XYOG)/nrowy
kroxx <- Reduce(`+`, kro)/nrowy
Bhat <- t(t(xyog)%*%kroxx)
BB <- ks::invvec(Bhat, ncol = (m*ncoly), nrow = (ncolylag + q))
resi <- list()
fitte <- matrix(nrow = nrowy, ncol = ncoly)
Ehat1 <- matrix(NA, ncol = ncoly, nrow = nrowy)
for (i in 1:nrowx){
resi[[i]] <- y[i, ] - t(Gtilde[[i]])%*%t(BB)%*%x[i,]
fitte[i,] <- t(t(Gtilde[[i]])%*%t(BB)%*%x[i,])
}
Ehat1 <- t(do.call("cbind", resi))
Omegahat <- (t(Ehat1)%*%Ehat1)/(nrowy)
data$Omegahat = Omegahat
data$BB = BB
#3. Convergence check
ll1 <- loglike(param = param1$par, data)
if((ll1-ll0)>0){
param <- as.matrix(matrixcalc::vec(cgam1))
}
err <- abs(ll1 - ll0)
ll0 <- ll1
loglik1[iter] <- ll1
message(paste("iteration", iter, "complete\n"))
cat(paste('Log-likelihood:', round(ll0,3), '\n'))
if (err<epsi | iter == n.iter) message('Converged\n')}
} else{
message('NLS estimation\n')
#Inizialization of iter
#Convergence algorithm
#1. NLS estimation of gamma and c with NLS estimates of Bhat and Omegahat
#2. NLS of Bhat with new values of gamma and c
#3. Convergence check
#4. 1-2-3 until convergence
while (iter < n.iter & err > epsi){
Sys.sleep(0)
iter <- iter+1L
#Parameters
low1 <- replicate(ny, 0)
#1.Maximum likelihood estimation of gamma and c
cl <- parallel::makeCluster(ncores) # set the number of processor cores
parallel::setDefaultCluster(cl=cl)
param1 <- optimParallel::optimParallel(par = as.vector(param), fn = SSQ, lower = c(low1, apply(y, 2, min)),
data = data, parallel = list(cl = cl, forward = FALSE, loginfo = FALSE))
parallel::stopCluster(cl)
cgam1 <- matrix(param1$par, ncol = 2L, nrow = (ny*(m-1)))
#2.NLS estimation of Bhat with new values of gamma and c
glog <- rep(0, ncol(y))
GT <- list()
Gtilde <- list()
kro <- list()
for (i in 1:nrowx){
for(t in 1:(m-1)){
for (j in 1 : ny){
glog[j] <- (1+exp(-cgam1[((t-1)*ny + j),1]*(st[i]-cgam1[((t-1)*ny + j),2])))^(-1)
}
if(singlecgamma == TRUE){
GT[[t]] <- diag(rep(glog[1], ncoly))
}else{
GT[[t]] <- diag(glog)
}
}
Gtilde[[i]] <- t(cbind(In, do.call(cbind,GT)))
kro[[i]] <- kronecker(Gtilde[[i]], x[i,])
}
Bhat <- MASS::ginv(t(t(do.call("cbind", kro)))%*%(t(do.call("cbind", kro))))%*%t(t(do.call("cbind", kro)))%*%(matrixcalc::vec(t(y)))
BB <- ks::invvec(Bhat, ncol = (m*ncoly), nrow = (ncolylag + q))
resi <- list()
resiresi <- list()
fitte <- matrix(nrow = nrowy, ncol = ncoly)
for (o in 1:nrowx){
resi[[o]] <- y[o, ] - t(Gtilde[[o]])%*%t(BB)%*%x[o,]
resiresi[[o]] <- resi[[o]]%*%t(resi[[o]])
fitte[o,] <- t(t(Gtilde[[o]])%*%t(BB)%*%x[o,])
}
Ehat1 <- Reduce("+", resiresi)
Omegahat <- Ehat1/(nrowy-1L)
data$Omegahat = Omegahat
data$BB = BB
#3. Convergence check
ll1 <- loglike(param = param1$par, data)
if((ll1-ll0)>0){
param <- as.matrix(matrixcalc::vec(cgam1))
}
err <- abs(ll1 - ll0)
ll0 <- ll1
loglik1[iter] <- ll1
message(paste("iteration", iter, "complete\n"))
cat(paste('Log-likelihood:',round(ll0, 3), '\n'))
if(iter > 50){
if(loglik1[[iter]] == loglik1[[(iter-2)]]){
iter <- n.iter
}}
if (err<epsi | iter == n.iter) message('Converged\n')}
}
#Calculating residuals and estimating standard errors
residuals1 <- t(do.call("cbind", resi))
varhat <- diag(Omegahat)
bb1 <- BB[,1:ncoly]
bb2 <- list()
for(t in 1:(m-1)){
bb2[[t]] <- as.data.frame(BB[,(ncoly*(t-1)+1):(ncoly*t)] +
BB[,(ncoly*(t)+1):(ncoly*(t+1))])
}
bb4 <- as.matrix(do.call(rbind,bb2))
BBhat <- rbind(bb1, bb4) ##Estimates of the coefficients
##Calculating standard errors, t-test and p-values
colnames(cgam1) <- c('gamma', 'c')
covbb <- matrix(nrow = m*ncolx, ncol = ncoly)
ttest <- matrix(nrow = m*ncolx, ncol = ncoly)
pval <- matrix(nrow = m*ncolx, ncol = ncoly)
ee <- matrix(nrow = m*ncolx, ncol = ncoly)
for (j in 1 : ncoly){
covbb[,j] <- sqrt(diag(MASS::ginv(t(x) %*%x))*varhat[j])
ttest[,j] <- BBhat[,j]/covbb[,j]
pval[,j] <- 2*pt(abs(ttest[,j]),df=(nrowy*m-m*(ncolx)), lower.tail = FALSE)
}
##Significances used for the summary function
signifi <- matrix('',nrow = nrow(pval), ncol = ncol(pval))
for(i in 1:nrow(pval)){
for(j in 1:ncol(pval)){
if(pval[i,j] < 0.01)
{signifi[i,j] <- '***'}
else if(pval[i,j] <=0.05 & pval[i,j]>0.01){
signifi[i,j] <- '**'
}
else if(pval[i,j] <= 0.1 & pval[i,j]>0.05){
signifi[i,j] <- '*'
}
}
}
##Calculating univariate log-likelihoods, AIC and BIC criteria
residui <- residuals1
omega1 <- diag(Omegahat)
k <- nrow(BBhat)
ll2 <- NULL
AIC1 <- NULL
BIC1 <- NULL
for (l in 1:ncoly){
ll2[l] <- -(nrow(y)/2)*log(omega1[l]) - (t(residui[,l])%*%residui[,l])/(2*omega1[l])
AIC1[l] <- 2*k - 2*ll2[l]
BIC1[l] <- -2*ll2[l] + k*log(nrowy)
}
names1 <- list()
for(j in 1:m){
names1[[j]] <- as.data.frame(paste(colnames(x), ' m_', j, sep = ''))
}
names1 <- as.matrix(do.call(rbind,names1))
rownames(BBhat) <- names1
colnames(BBhat) <- colnames(y)
modeldata <- list(y, x)
fitte <- fitte[!is.na(fitte[,1]),]
results <- list(BBhat, covbb, ttest, pval, cgam1, Omegahat, fitte, residuals1, ll1, ll2, AIC1, BIC1, Gtilde, modeldata, BB, m, p,
st, y, exo, constant, method, singlecgamma)
names(results) <- c('Bhat','StDev', 'ttest', 'pval', 'Gammac', 'Omega', 'fitted', 'residuals', 'MultiLL', 'LL', 'AIC',
'BIC', 'Gtilde', 'Data', 'B', 'm', 'p', 'st', 'yoriginal', 'exo', 'constant', 'method', 'singlecgamma')
class(results) = 'VLSTAR'
return(results)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.