#VAR(p) implementation based on Bulteel et al. 2018, modified by Tim Leers
#shiny-specific arguments needed for computedata
computeDataArgs.var <- function(model){
list(input$nVar,
input$nTime,
0,
input$selection1,
val=TRUE,
burn=1000,
#model-specific parameters
current_phi_input(),
current_inno_input())
}
modelDataArgs.var <- function(model){
return(
list(
input$selection1,
filedata_updated(),
selectedLagNum(),
loaded_dataset_index_variable())
)
}
relevantModelParameters.varest<-function(tmod){
return(list(phi=extractPhi(tmod),
inno=extractInno(tmod)
))
}
currentModelParameters.varest <-function(tmod){
return(list(phi=current_phi_input(),
inno=current_inno_input())
)
}
currentModelParameters.var <- currentModelParameters.varest
relevantModelParameters.var <- relevantModelParameters.varest
####Name function-----
modelName.var<-function(model){
return('Vector Autoregression')
}
modelName.varest<-modelName.var
####Data-generating function----
#SOURCE NOTE: Code is based on code distributed by the paper of Bulteel, Tuerlinckx, Brose, and Ceulemans
#TITLE:Improved Insight into and Prediction of Network Dynamics by Combining VAR and Dimension Reduction
#doi:10.1080/00273171.2018.1516540
computeData.var <-function(nVar,
time,
error,
model,
val=TRUE,
burn=1000,
mod_vars,
...){
inno <- mod_vars$inno
phi <- mod_vars$phi
#To avoid crashes, we validate the phi matrix and the innovation matrix.
if(val){
if(!validate_phi(phi)){
warning("Transition matrix invalid")
return(NULL)
}
if(!validate_inno(inno)){
warning("Innovation matrix invalid")
return(NULL)
}
print("Matrices are valid.")
}
#Generate errors
innovations <- mvtnorm::rmvnorm(time+burn,rep(0,nVar),inno)
#Create empty matrix
U <- matrix(innovations, time + burn, nVar)
simdata <- matrix(0, time + burn, nVar)
simdata[1, ] <- U[1, ]
#withProgress(message = paste0('Simulating ',model), value = 0, {
for (row in 2:(time + burn)) {
simdata[row, ] = simdata[(row - 1), ] %*% phi + U[row, ]
}
randomError <- matrix(rnorm(time * nVar, 0, 1), time, nVar)
#E <- sqrt(error) * randomError
E <- 0
Y <- simdata[-(1:burn), ] + E
#})
return(Y)
}
computeData.varest<-computeData.var
output$var_sim_output <- renderUI({
tagList(
transitionMatrixUI(ns(session)$ns,'phi'),
innovationMatrixUI(ns(session)$ns,'inno')
)
})
output$var_mod_output <- renderUI({
})
prev_mod <<- reactiveVal(NULL)
prev_mod2 <<- reactiveVal(NULL)
modelDataParams.var<-function(model){
}
####Model fit function----
modelData.varest <- function(model, dataset, lagNum, index_vars = NULL,...) {
require(vars)
# variables<-na.omit(str_match(names(dataset),id))
# day_index <- 'UNIT'
# beep_index <- 'OCCASION'
# mlVAR(dataset,variables,id,lags=1,
# dayvar=day_index,
# beepvar=beep_index,
# estimator='default',
# temporal='correlated')
if(!is.null(index_vars)) {
VAR(dataset, type = "const", p = lagNum, exogen = dataset[,index_vars])
} else {
VAR(dataset, type = "const", p = lagNum)
}
# ar(
# dataset,
# type = 'const',
# p = lagNum,
# method = 'burg',
# aic=FALSE,
# order.max=ncol(dataset)*ncol(dataset)
#
}
modelData.var <- modelData.varest
####Parameter extraction functions-----
extractPhi.varest <- function(model){
nvar <- ncol(model$y)
VAR_coeff <- matrix(0, nvar, nvar)
for (i in 1:length(model$varresult)) {
VAR_coeff[i, ] <-
model$varresult[[i]]$coefficients[1:nvar]
}
phi <- VAR_coeff
}
extractInno.varest <- function(model){
nvar <- ncol(model$y)
resid_matrix <- matrix(0,length(model$varresult[[1]]$residuals),nvar)
VAR_coeff <- matrix(0, nvar, nvar)
for (i in 1:length(model$varresult)) {
resid_matrix[,i] <- model$varresult[[i]]$residuals
}
return(cov(na.omit(resid_matrix)))
}
####Residual extraction function----
extractResiduals.var <- function(model){
residuals(model)
}
extractResiduals.varest <- extractResiduals.var
####Error computation function----
computeError.varest <- function(model,pred,dat){
error <- pracma::rmserr(dat %>% unlist() %>% as.numeric(), pred[[1]] %>% as.numeric())
sqr_resid <- (dat-pred[[1]])^2
#sd <- apply(sqr_resid,2,sd)/sqrt(nrow(dat))
sd <- sd(as.numeric(unlist(sqr_resid)))/sqrt(nrow(dat))
error<-rlist::list.append(error,sd)
names(error)[[7]] <- 'sd'
return(error)
}
computeError.var <- computeError.varest
####Model prediction function----
altpredict.varest <- function(model,data){
pred <- matrix(0,nrow(data),ncol(data))
#se <- matrix(0,n.ahead,ncol(data))
se <- NULL
data <- matrix(unlist(data),ncol=ncol(extractResiduals.var(model)))
for (i in 1:ncol(data)){
# tmp<-try(suppressWarnings(stats:::predict.ar(model$fit_list[[i]],data[,i],n.ahead=n.ahead)))
# #print(class(tmp))
# if(class(tmp)!='list'){
# tmp<-stats:::predict.ar(model$fit_list[[i]],data[,i],n.ahead=n.ahead,method='burg')
# warning('predict.ar failed using ols - using burg instead')
# }
#withProgress(message = 'Simulating data', value = 0, {
phi <- extractPhi(model)
pred[1,] <- data[1,]
for (row in 2:nrow(data)) {
pred[row, ] = phi %*% data[row-1,]
}
#se[,i] <- as.vector(tmp$se)
}
return(list(pred,se))
}
altpredict.var<-altpredict.varest
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.