[R] : Extra loops in order to fill in a 3D table

[Harris Agiropoulos]

Hello all,

I would like to have two or maybe three extra “for”. One for changing beta_x with values 0.00 and 0.20 and the others for changing phi_x and phi_y with values 0.00, 0.30, 0.90.
Does anyone know how to implement this?

library(MASS)
library(forecast) 
library(lmtest)
library("dyn")
library(nlme)
mysummary <- function(betax,npar=TRUE,print=TRUE) {
  options(scipen=999)
  set.seed(1234)
  par(mfrow=c(1,2))
  rep<-10
  n<-100
  sigma<-1
  mx<-0.8
  my<-0.8
  theta_y<-2
  theta_x<-2
  beta_x<-0.00
  beta_y<-betax
  phi_x<-0.00
  phi_y<-0.00
  break1<-n/2
  break2<-n/2
  break3<-n/5
  break4<-n/5
  vhta<-c(NA,rep)
  alpha<-c(NA,rep)
  da.wa<-c(NA,rep)
  tau0<-c(NA,rep)
  tau<-c(NA,rep)
  p_val<-c(NA,rep)
  vhta2<-c(NA,rep)
  tau2<-c(NA,rep)
  tau3<-c(NA,rep)
  tau4<-c(NA,rep)
  tau5<-c(NA,rep)
  tau6<-c(NA,rep)
  p_val2<-c(NA,rep)
  time<-c(1:n)
  for (i in 1:rep)
  {
    t <- c(1:n) # specifying the time
    dummy <- as.numeric(ifelse(t <= break1, 0, 1)) # specifying the dummy for trend break at T = 40
    dummy2 <- as.numeric(ifelse(t <= break2, 0, 1)) # specifying the dummy for trend break at T = 80
    dummy3 <- as.numeric(ifelse(t <= break3, 0, 1)) # specifying the dummy for trend break at T = 20
    dummy4 <- as.numeric(ifelse(t <= break4, 0, 1)) # specifying the dummy for trend break at T = 70
    z<-w<-rnorm(n,sd=sigma)
    for (t in 2:n) z[t]<-phi_y*z[t-1]+w[t]
    Time<-1:n
    ar1_1 <- ts(my + beta_y*Time - theta_y*dummy + z) #- theta_y*(Time-dummy2)*dummy2 + z) # This is the trend stationary model with break in trend
    y <- ts(my + beta_y*Time - theta_y*dummy) #- theta_y*(Time-dummy2)*dummy2) # This is just the trend line that we see in "red" in the plot below
    plot(ar1_1, main = "Simulated series with breaks at T = 40,80")
    lines(y, col = "red") ## Plotting a sample of the model that we have simulated
    ## Now we will simulate the sample data above 1000 times and check for unit roots for each of these samples ##
    g<-s<-rnorm(n,sd=sigma)
    for (t in 2:n) g[t]<-phi_x*g[t-1]+s[t]
    Time<-1:n
    ar1_2 <- ts(mx + beta_x*Time - theta_x*dummy3 +g) #- theta_x*(Time-dummy4)*dummy4 + g) # This is the trend stationary model with break in trend
    x <- ts(mx + beta_x*Time - theta_x*dummy3) #- theta_x*(Time-dummy4)*dummy4) # This is just the trend line that we see in "red" in the plot below
    plot(ar1_2, main = "Simulated series with breaks at T = 20,70")
    lines(x, col = "red") # Plotting a sample of the model that we have simulated
    lmold.r<-lm(ar1_1 ~ ar1_2)
    lm.r<-lm(ar1_1 ~ Time + ar1_2)
    Data<-data.frame(ar1_1, ar1_2, Time)
    fglm.r<-gls(ar1_1 ~ Time + ar1_2, Data, corAR1(0.9,form = ~ 1,TRUE))
    fglm2.r<-gls(ar1_1 ~ Time + ar1_2, Data, corAR1(0.5,form = ~ 1,TRUE))
    fglm3.r<-gls(ar1_1 ~ Time + ar1_2, Data, corAR1(0.3,form = ~ 1,TRUE))
    fglm4.r<-gls(ar1_1 ~ Time + ar1_2, Data, corAR1())
    fglm5.r<-gls(ar1_1 ~ Time + ar1_2, Data, weights = varIdent(form =~ 1))
    tau0[i]<-coef(summary(lmold.r))["ar1_2","t value"] #simple linear regression 
    tau[i]<-coef(summary(lm.r))["ar1_2","t value"] #simple linear regression with trend
    tau2[i]<-summary(fglm.r)$tTable["ar1_2","t-value"] #GLS with rho = 0.9
    tau3[i]<-summary(fglm2.r)$tTable["ar1_2","t-value"] #GLS with rho = 0.5
    tau4[i]<-summary(fglm3.r)$tTable["ar1_2","t-value"] #GLS with rho = 0.3
    tau5[i]<-summary(fglm4.r)$tTable["ar1_2","t-value"] #GLS with rho_hut
    tau6[i]<-summary(fglm5.r)$tTable["ar1_2","t-value"] #GLS with variance power
  }
  tau0
  tau
  tau2
  tau3
  tau4
  tau5
  tau6
  ti0<-sum(abs(tau0) > 1.96)
  ti<-sum(abs(tau) > 1.96)
  ti2<-sum(abs(tau2) > 1.96)
  ti3<-sum(abs(tau3) > 1.96)
  ti4<-sum(abs(tau4) > 1.96)
  ti5<-sum(abs(tau5) > 1.96)
  ti6<-sum(abs(tau6) > 1.96)
  results<-cbind(c(t0=ti0,t=ti,t5=ti5,t4=ti4,t3=ti3,t2=ti2,t6=ti6))
  results
  
}

betax<-c(0.00,0.20)
x<-sapply(betax,mysummary)
x


	[[alternative HTML version deleted]]