R/impulse_response.R

In vwendy/pompom: Person-Oriented Method and Perturbation on the Model

# 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),]
  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)
}