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]]
______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
