Nothing
# generate impulse response (multivariate time-series) from beta matrix
# return value in long format, each column is where the impulse goes into,
#last column is number of time steps
# so the data organization return value is:
# long format
# var.number by var.number columns
# the 1st var.number columns are the impulse reseponse received by the 1st variables
# the 2nd var.number columns are the impulse reseponse received by the 2nd variables
# ...
impulse.response <- function(var.number,
lag.order = 1,
steps = 100,
beta.matrix)
{
beta <- beta.matrix
n.obs <- steps
# convert to A matrix and Phi matrix
alpha <- as.matrix(beta[(var.number+1):(2*var.number), (var.number+1):(2*var.number)],
nrow = var.number,
ncol = var.number,
byrow = TRUE)
phi <- as.matrix(beta[(var.number+1):(2*var.number), 1:var.number],
nrow = var.number,
ncol = var.number,
byrow = TRUE)
identity.matrix <- diag(1, var.number, var.number)
# by algebra, coefficient matrix is the inverse of (I -A) * Phi
coeff.matrix.eta.lag <- solve(identity.matrix-alpha) %*% phi
return.time.series.shock <- NULL
return.time.series.shock$steps <- 1:n.obs
for (node.to.shock in 1:var.number)
{
shock <- rep(0,var.number)
shock[node.to.shock] <- 1
shock <- as.matrix(shock, nrow = var.number, ncol = 1)
time.series.shock <- matrix(rep(0, var.number * steps), nrow = var.number, ncol = steps)
start.point <- 1 # at least 2, to give timepoint 1 as initial value
for(col in 1:start.point)
{
time.series.shock[,col] <- solve(identity.matrix-alpha) %*% shock
}
# noise come as psi matrix directs
for (time in (start.point+1):n.obs)
{
# ONLY shock at T =1
time.series.shock[,time] <- coeff.matrix.eta.lag %*%
time.series.shock[,time-1]
}
column.names <- paste("from.", node.to.shock, ".to.", 1:var.number, sep = "")
time.series.shock <- data.frame(t(time.series.shock))
names(time.series.shock) <- column.names
return.time.series.shock <- cbind(return.time.series.shock,
time.series.shock)
}
return(return.time.series.shock)
}
# compute recovery time from a time.series
compute.recovery.time.ts <- function(time.series, steps, var.number, threshold = 0.01){
recovery.time <- 0
for (row in steps:1)
{
if ((time.series[row] >= threshold ||
time.series[row] <= -threshold)
&& recovery.time ==0)
{
recovery.time <- row + 1
}
}
return (recovery.time)
}
compute.recovery.time.from.beta <- function(beta.matrix,
var.number,
lag.order =1,
threshold = 0.01,
replication = 200,
steps = 100){
recovery.time <- matrix(rep(0, var.number * var.number), nrow = var.number)
time.profiles <- impulse.response(var.number,
lag.order,
steps,
beta.matrix)
index <- 2 # first column is the time steps
for (from in 1:var.number)
{
for (to in 1:var.number)
{
recovery.time[from,to] <- compute.recovery.time.ts(time.profiles[,index], steps, var.number, threshold)
index <- index + 1
}
}
return (recovery.time)
}
# To generate beta matrices based on point estimation, SE
# each row is one beta estimation link (according to the output lavaan)
# first 4 columns are edge to (lhs from lavvan), edge from (rhs from lavaan), estimation of beta, SE of beta
# after the first 4, each column is a replication generated by the estimation and SE
bootstrap.beta.from.beta.est <- function(beta.est, replication = 200, var.number){
# rm(.Random.seed, envir=.GlobalEnv)
set.seed(1234)
bootstrapped.beta <- matrix(rep(0, (replication + 4) * length(beta.est[,1])),
nrow = length(beta.est[,1]),
ncol = replication + 4,
byrow = TRUE)
# replication <- 1000 # debug
for (row in 1:length(beta.est[,1]))
{
beta.bootstrap <- rnorm(replication, mean = beta.est$est[row], sd = beta.est$se[row])
bootstrapped.beta[row,1] <- beta.est$lhs[row]
bootstrapped.beta[row,2] <- beta.est$rhs[row]
bootstrapped.beta[row,3] <- beta.est$est[row]
bootstrapped.beta[row,4] <- beta.est$se[row]
bootstrapped.beta[row,5:(replication + 4)] <- beta.bootstrap
}
return(bootstrapped.beta)
}
# to parse one replication (repnum) in the returned value from "bootstrap.beta.from.list" into a beta matrix
# size of matrix is nrow = (lag.order+1) * var.number, ncol = nrow
parse.bootstrapped.beta.to.beta.matrix <- function(var.number,
replication,
boostraped.beta,
repnum,
lag.order) # pass the modelfit (lavaan)
{
beta.est <- data.frame(cbind(boostraped.beta[,1:2],boostraped.beta[,4 + repnum]))
names(beta.est) <- c("lhs", "rhs", "est")
beta.est$est <- as.numeric(as.character(beta.est$est))
beta.matrix <- matrix(rep(0, (var.number*(lag.order + 1))^2),
nrow = (var.number*(lag.order+1)),
ncol = (var.number*(lag.order+1)),
byrow = TRUE)
lhs.number <- as.numeric(gsub("eta","", beta.est$lhs))
rhs.number <- as.numeric(gsub("eta","", beta.est$rhs))
for(i in 1:length(beta.est[,1]))
{
beta.matrix[lhs.number[i], rhs.number[i]] <- beta.est$est[i]
}
return (beta.matrix)
}
# function to simulate time profiles (impulse response analysis) by beta (bootstrapped)
# the data organization is each draw of beta has a full var.number by var.number time profiles, from function impulse.response
# for each impulse.response output the data organization return value is:
# long format
# var.number by var.number columns
# the 1st var.number columns are the impulse reseponse received by the 1st variables
# the 2nd var.number columns are the impulse reseponse received by the 2nd variables
# ...
# then this is attached with a column of replication number
# and rbind with other replications
# then each column has N = (time profile steps) * (replication number) rows
# each column is also a node-to-node impulse response with multiple replications
# you can then compute the mean the confidence interval of the impulse response by columns
impulse.response.boot <- function(replication = 200,
var.number,
lag.order,
bootstrapped.beta,
steps)
{
simulation <- NULL
for (repnum in 1:replication)
{
one.beta.draw <- parse.bootstrapped.beta.to.beta.matrix(var.number,
replication,
bootstrapped.beta,
repnum,
lag.order)
impulse.response.time.series <- impulse.response(var.number,
lag.order,
steps,
one.beta.draw)
impulse.response.time.series$repnum <- rep(repnum,
nrow(impulse.response.time.series))
simulation <- rbind(simulation, impulse.response.time.series)
}
return (simulation)
}
compute.confidence.interval <- function(dataset, mean)
{
n.obs <- length(dataset)
cutoff.obs <- round(n.obs * .025,0)
dataset <- data.frame(dataset)
# dataset <- dataset[order(dataset),] # code before 1/22/2021
# print("inside corrected code!")
# print(str(dataset))
names(dataset) <- "recovery.time" # code after 1/22/2021
dataset <- dataset[order(dataset$recovery.time),] # code after 1/22/2021
lower <- dataset[cutoff.obs + 1]
upper <- dataset[n.obs - cutoff.obs + 1]
confidence.interval <- list(lower, upper)
return(list(lower = lower, upper = upper))
}
iRAM.boot <- function(model.fit,
var.number,
lag.order = 1,
threshold = 0.01,
replication = 200,
steps = 100
)
{
recovery.time.reps <- matrix(rep(0, var.number*var.number*replication),
nrow = replication,
ncol = var.number*var.number)
recovery.time.mean <- matrix(rep(0, var.number*var.number),
nrow = 1,
ncol = var.number*var.number)
recovery.time.ci.upper <- matrix(rep(0, var.number*var.number),
nrow = 1,
ncol = var.number*var.number)
recovery.time.ci.lower <- matrix(rep(0, var.number*var.number),
nrow = 1,
ncol = var.number*var.number)
ret.mean <- matrix(rep(0, var.number*var.number), nrow = var.number)
ret.upper <- matrix(rep(0, var.number*var.number), nrow = var.number)
ret.lower <- matrix(rep(0, var.number*var.number), nrow = var.number)
# parse out beta estimates
beta.est <- parse.beta.as.estimates(var.number, model.fit)
beta.reps <- bootstrap.beta.from.beta.est(beta.est, replication, var.number)
# simulate by model fit
simulation.data <- impulse.response.boot(replication, var.number, lag.order, beta.reps, steps)
# compute recovery time per column and replication
for (from.node in 1:var.number)
{
for (to.node in 1:var.number)
{
# compute the column number in the simulated dataset
colnum <- (from.node-1)*var.number + to.node + 1 # there is an extra steps column in simulation data
for (repnum in 1:replication)
{
row.in.simu <- (repnum -1) * steps + 1
recovery.time.reps[repnum, colnum -1] <- compute.recovery.time.ts(
simulation.data[row.in.simu:(row.in.simu+steps-1), colnum], steps, var.number, threshold)
}
}
}
# compute mean and confidence interval by each node-to-node cell
recovery.time.mean <- colMeans(recovery.time.reps)
for (col in 1:(var.number*var.number))
{
recovery.time.ci.lower[1,col] <- compute.confidence.interval(recovery.time.reps[,col])$lower
recovery.time.ci.upper[1,col] <- compute.confidence.interval(recovery.time.reps[,col])$upper
}
# adjust to matrix form
for (from in 1:var.number)
{
for(to in 1:var.number)
{
ret.mean[from, to] <- recovery.time.mean[(from-1)*var.number + to]
ret.lower[from, to] <- recovery.time.ci.lower[1,(from-1)*var.number + to]
ret.upper[from, to] <- recovery.time.ci.upper[1,(from-1)*var.number + to]
}
}
# return value should be a mean and confidence interval vector with length of var.number by var.number
return(list(mean = ret.mean,
lower = ret.lower,
upper = ret.upper,
time.profile.data = simulation.data,
recovery.time.reps = recovery.time.reps))
}
#' Generate iRAM (impulse response anlaysis metric) from model fit.
#'
#'
#' @param model.fit model fit object generated by lavaan
#' @param beta beta matrix for a point estimate
#' @param var.number number of variables in the time series
#' @param lag.order lag order of the model to be fit
#' @param threshold threshold of calculation of recovery time (duration of perturbation), default value is 0.01
#' @param boot to bootstrap, default value is FALSE
#' @param replication number of replication of bootstrap, default value is 200
#' @param steps number of steps of impulse response analysis, default value is 100
#'
#' @return iRAM matrix. Rows represent where the orthognal impulse was given, and columns represent the response. Dimension is var.number by var.number.
#'
#'
#' @references Lütkepohl, H. (2007). New introduction to multiple time-series analysis. Berlin: Springer.
#'
#' @examples
#' \dontshow{
#' boot.iRAM <- iRAM(model.fit = usemmodelfit,
#' beta = NULL,
#' var.number = 3,
#' lag.order = 1,
#' threshold = 0.01,
#' boot = TRUE,
#' replication = 50, # default replication is 200, reduced to 50 to shorten running time
#' steps = 30 # default steps is 100, reduced to 30 to shorten running time
#' )
#'boot.iRAM$mean
#' }
#' \donttest{
#' boot.iRAM <- iRAM(model.fit = usemmodelfit,
#' beta = NULL,
#' var.number = 3,
#' lag.order = 1,
#' threshold = 0.01,
#' boot = TRUE,
#' replication = 200,
#' steps = 100
#' )
#'boot.iRAM$mean
#' }
#'
#'
#' @export
#'
iRAM <- function(model.fit,
beta,
var.number,
lag.order = 1,
threshold = 0.01,
boot = FALSE,
replication = 200,
steps = 100)
{
if (!is.null(model.fit))
{
if (boot)
{
iRAM.boot.rt <- iRAM.boot(model.fit = model.fit,
var.number = var.number,
lag.order = lag.order,
threshold = threshold,
replication = replication,
steps = steps
)
return (iRAM.boot.rt)
}
else
{
beta.matrix <- parse_beta(var.number,
model.fit,
lag.order,
matrix = T)
time.series.data <- impulse.response(var.number,
lag.order,
steps = steps,
beta.matrix$est)
# print(beta.matrix)
rt <- compute.recovery.time.from.beta(beta.matrix$est,
var.number,
lag.order,
threshold ,
replication,
steps)
return(list(recovery.time = rt,
time.series.data = time.series.data))
}
}
else
{
# this is when you only have true beta matrix, no model fit
# you can use the beta matrix (second parameter) to get the true value of iRAM
# rt <- compute.recovery.time.from.beta(beta,
# var.number,
# lag.order,
# threshold ,
# replication,
# steps)
# return(list(recovery.time = rt))
# beta.matrix <- beta
time.series.data <- impulse.response(var.number,
lag.order,
steps = steps,
beta)
# print(beta.matrix)
rt <- compute.recovery.time.from.beta(beta,
var.number,
lag.order,
threshold ,
replication,
steps)
return(list(recovery.time = rt,
time.series.data = time.series.data))
}
}
#' Generate iRAM (impulse response anlaysis metric) in the equilibrium form.
#'
#'
#' @param beta.matrix beta matrix for a point estimate
#' @param var.number number of variables in the time series
#' @param lag.order lag order of the model to be fit
#'
#' @return a list of equilibria. First numeric number in the variable name indicate where the impulse was given, and the second numeric number indicate the response, e.g., e12 indicates equilibrium of node 2 when node 1 is given an impulse.
#'
#' @examples
#' \dontshow{
#' iRAM_evalue <- iRAM_equilibrium(beta.matrix = true_beta_3node,
#' var.number = 3,
#' lag.order = 1
#' )
#' iRAM_evalue
#' }
#' \donttest{
#' iRAM_evalue <- iRAM_equilibrium(beta.matrix = true_beta_3node,
#' var.number = 3,
#' lag.order = 1
#' )
#' iRAM_evalue
#' }
#'
#'
#' @export
#'
#'
iRAM_equilibrium <- function(beta.matrix, var.number, lag.order)
{
phi <- as.matrix(beta.matrix[(var.number+1):(2*var.number),1:var.number]) # lag-1 relations
alpha <- as.matrix(beta.matrix[(var.number+1):(2*var.number),(var.number+1):(2*var.number)]) # contemporaneous relations
coeff.mat <- solve(diag(var.number)-alpha) %*% phi
equilibrium.mat <- coeff.mat %*% solve(diag(var.number) - coeff.mat) %*%
solve(diag(var.number)-alpha) %*% diag(var.number)
equilibrium.list <- t(as.vector(equilibrium.mat))
equilibrium.list <- data.frame(equilibrium.list)
for (first.index in 1:var.number){
for(second.index in 1:var.number)
{
names(equilibrium.list)[ var.number * (first.index - 1) + second.index] <-
paste("e", first.index, second.index, sep = "")
}
}
return (equilibrium.list)
}
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.