Here is my trended R code to plot trended vs detrended data

test.R

## Stephen J. Barr
##   - this is some practice retrending the model

setwd("/home/stevejb/myhg/is-solver/indep_sim/")


retrendDemo <- function(THETA = 0.7880) {

  ## produce initial values
  x = rnorm(30, sd=.03)
  K = 500

  z0 = 1
  zmat = matrix(0, ncol = length(x), nrow=1)
  zmat[1,1] = z0

  theta = 0.788

  ztrend = (z0 * cumprod(exp(x)))^(1/(1-theta))

  times = 1:length(x)

  ########################################
  ## Plot comparaing the trended (blue)
  ## vs untrended (red)
  plot(times, ztrend, type="l", col="blue",
       ylim=(c(-2, max(ztrend))),
       main = "Trended vs Untrended Shock")

  par(new=T)

  plot(times, x, type="l", col="red",
       xaxt='n', axes=F)
  axis(4, pretty(c(min(x), max(x))), col='red')
  grid()
  ########################################

  dev.new()

  Kuntrended = matrix(K, ncol = length(x), nrow=1)
  Ktrended = matrix(K, ncol = length(x), nrow=1)

  ## K_t = k_t * {z^P_{t-1}}^(1/(1-THETA))
  for(t in 2:length(x)) {
    Ktrended[1,t] = K * ztrend[(t-1)]^(1/(1-THETA))
  }


  meaninv = mean((Ktrended[2:30] - 0.85*(Ktrended[1:29]))/Ktrended[1:29])
  plot(times, log(Ktrended), type="l",
       main = paste("Mean inv: ", meaninv))

  print("Press enter to continue")
  y <- scan(n=1)
  dev.off()
  dev.off()
  
}


retrendFromData <- function(
                     EPS_P_PATH = "bv9___--is--tr--epsilon_P_RAW.csv",
                     K_PATH = "bv9___--is--tr--k_prev_raw.csv",
                     P_PATH = "bv9___--is--tr--p_prev_raw.csv",
                     doplot = FALSE
                     ) {
  

  ### VARS
  THETA = 0.7880
  DELTA = 0.15
  z0 = 1
  #########


  epsP = read.csv(EPS_P_PATH, header=FALSE, as.is=TRUE)
  kmat = read.csv(K_PATH, header=FALSE, as.is=TRUE)
  pmat = read.csv(P_PATH, header=FALSE, as.is=TRUE)

  epsP = as.matrix(epsP)
  kmat = as.matrix(kmat)
  pmat = as.matrix(pmat)

  ########################################
  ## trended k plot
  if(doplot) {
    mycolors = rainbow(10)
    for(t in 1:min(nrow(kmat),10)) {
      if(t == 1) {
        plot(kmat[t,], type="l", col = mycolors[t],
             main="DETRENDED K EXAMPLES")
      } else {
        lines(kmat[t,], type="l", col = mycolors[t])
      }
    } ## end for
  }
  ########################################


  
  ## STEP 1 - CALCULATE THE CUMULATIVE PERMANENT SHOCK
  ## STEP 2 - RETREND THE K MATRIX
  ## STEP 3 - RETREND THE P MATRIX

  ## STEP 1 -- 
  ztrend_mat = matrix(0, nrow=nrow(epsP),
    ncol=ncol(epsP))  
  for( nn in 1:nrow(epsP)) {
    ztrend_mat[nn,] = cumprod(exp(epsP[nn,]))##^(1/(1-THETA))
  }

  ## STEP 2 -- RETREND INV
  Ktrended = matrix(0, ncol = ncol(kmat),
    nrow = nrow(kmat))
  Ktrended[,1] = kmat[,1]
  for(nn in 1:nrow(kmat)) {
    for(tt in 2:ncol(kmat)) {
      Ktrended[nn,tt] = kmat[nn,tt] * (ztrend_mat[nn, (t-1)])^(1/(1-THETA))       
    }
  }

  ## STEP 3 -- RETREND DEBT
  Ptrended = matrix(0, ncol = ncol(pmat), nrow = nrow(pmat))    
  Ptrended[,1] = pmat[,1]
  for(nn in 1:nrow(kmat)) {
    for(tt in 2:ncol(kmat)) {
      Ptrended[nn, tt] = pmat[nn,tt] * ztrend_mat[nn,tt]
    }
  }
  d_indic_trended  = pmax(Ptrended, 0)
  leverage_trended = d_indic_trended/Ktrended

  d_indic_dtr = pmax(pmat, 0)
  leverage_dtr = d_indic_dtr/kmat
 
  make_inv_vector = function(kvec) {
    KL = length(kvec)
    invs = (kvec[2:KL] - (1-DELTA) * kvec[1:(KL-1)])/kvec[1:(KL-1)]
    return(invs)
  }

  ## CALCULATED TRENDED INVESTMENT
  ## 
  invmat_TR = t(apply(Ktrended, 1, make_inv_vector))
  invmat_DTR = t(apply(kmat, 1, make_inv_vector))  
  mean_inv_TR = mean(invmat_TR)
  mean_inv_DTR = mean(invmat_DTR)
  var_inv_TR = var(as.vector(invmat_TR))
  var_inv_DTR = var(as.vector(invmat_DTR))


  ## CALCULATE TRENDED LEVERAGE
  ##
  mean_lev_TR = mean(leverage_trended)
  mean_lev_DTR = mean(leverage_dtr)
  var_lev_TR = var(as.vector(leverage_trended))
  var_lev_DTR = var(as.vector(leverage_dtr))
  

  ## reslist = list(mean_TR = mean_TR,
  ##   mean_DTR = mean_DTR,
  ##   var_TR = var_TR,
  ##   var_DTR = var_DTR)
  ## return(reslist)
  print(paste("Mean inv TR", mean_inv_TR))
  print(paste("Mean inv DTR", mean_inv_DTR))
  print(paste("Mean lev TR", mean_lev_TR))
  print(paste("Mean lev DTR", mean_lev_DTR))
}

retrendFromData()