Hello,
You should post to the list, not just send your questions to me. If you
Cc the list the odds of getting more and better answers are bigger.
As for your question,
1. You forgot a variable X1, like this your code doesn't run. In what
follows I assume it's also a runif(40, 0, 20).
2. Some simplifications in your code could be made. I've left the
original instructions, commented out, and the corrections below them.
3. When you say that the bias and the se's are very small what do you
mean exactly? What kind of values are you expecting and what are the
results you're getting?
fun <- function(A) {
model <- lm(Y ~ X1 + X2, data = A)
e <- resid(model)
mylag <- function(e, d = 1) {
n <- length(e)
c(rep(NA, d), e)[1:n]
}
n <- length(e)
e1 <- mylag(e)
modele <- lm(e ~ e1 - 1)
rho <- coef(modele)
reps <- 20
for (i in 1:reps) {
n <- length(e)
x1star <- c(X1[1] * (1 - rho), X1[2:n] - rho * X1[1:(n - 1)])
x2star <- c(X2[1] * (1 - rho), X2[2:n] - rho * X2[1:(n - 1)])
ystar <- c(Y[1] * (1 - rho), Y[2:n] - rho * Y[1:(n - 1)])
#cr <- (1 - rho)
#cr[1:n] <- cr
cr <- rep(1- rho, n)
#W <- 1
#W[1:n] <- W
W <- rep(1, n)
W[1] = (1 - rho^2)/((1 - rho)^2)
#W <- c(W[1], W[2:n])
Model <- lm(ystar ~ cr + x1star + x2star - 1, weights = W)
bstar <- coef(Model)
#B0 <- (bstar[[1]][[1]])
#B1 <- bstar[[2]][[1]]
#B2 <- bstar[[3]][[1]]
#Yhat <- B0 + B1 * X1 + B2 * X2
#u <- Y - Yhat
u <- resid(Model)
myu <- function(u, d = 1) {
n <- length(u)
c(rep(NA, d), u)[1:n]
}
u1 <- myu(u)
modelu <- lm(u ~ u1 - 1)
Rho <- coef(modelu)
if (abs(Rho) > 1)
break
diff <- abs(Rho - rho)
if (diff < 1e-05)
break else rho <- Rho
}
return(coef(Model))
}
l <- runif(40, 0, 20)
X1 <- X2 <- runif(40, 0, 20)
U <- rnorm(1, 0, 2)
for (k in 2:40) {
U[k] = 0.9 * U[k - 1] + rnorm(1, 0, 1)
}
Y <- l + 2 * X1 + 5 * X2 + U
A <- data.frame(X1 = X1, X2 = X2, Y = Y)
result <- boot::tsboot(A, statistic = fun, R = 1000, l = nrow(A), sim =
"geom", orig.t = TRUE)
result
Call:
boot::tsboot(tseries = A, statistic = fun, R = 1000, l = nrow(A),
sim = "geom", orig.t = TRUE)
Bootstrap Statistics :
original bias std. error
t1* 12.563106 0.12260192 0.34212328
t2* 6.980546 -0.01316851 0.03659207
WARNING: All values of t3* are NA
Are these results ok?
Rui Barradas
Em 29-11-2012 04:34, Hock Ann Lim escreveu:
Dear Dr. Rui Barradas,
I have this following R programming code to stationary bootstrap the Cochrane-Orcutt Prais Winsten iterative method to obtain the parameter estimates.
I think there must be some problems in my programming as the bias and the std. error shown in the output are very small. Due to my little knowledge in R, I'm not able to rectify the problems. Hope you can help to point out the mistakes for me.
X1<-runif(40,0,20)
X2<-runif(40,0,20)
U<-rnorm(1,0,2)
for (k in 2:40){
U[k]=0.9*U[k-1]+rnorm(1,0,1)}
Y<-1+2*X1+5*X2+U
A <- data.frame(X1 = X1,X2=X2, Y = Y)
fun <- function(A){
model<-lm(Y~X1+X2,data=A)
e<-resid(model)
mylag<-function(e,d=1) {
n<-length(e)
c(rep(NA,d),e)[1:n]
}
n<-length(e)
e1<-mylag(e)
modele<-lm(e~e1-1)
rho<-coef(modele)
reps<-20
for(i in 1:reps){
n<-length(e)
x1star<-c(X1[1]*(1-rho),X1[2:n]-rho*X1[1:(n-1)])
x2star<-c(X2[1]*(1-rho),X2[2:n]-rho*X2[1:(n-1)])
ystar<-c(Y[1]*(1-rho),Y[2:n]-rho*Y[1:(n-1)])
cr<-(1-rho)
cr[1:n]<-cr
W<-1
W[1:n]<-W
W[1]=(1-rho^2)/((1-rho)^2)
W<-c(W[1],W[2:n])
Model<-lm(ystar~cr+x1star+x2star-1,weights=W)
bstar<-coef(Model)
B0<-(bstar[[1]][[1]])
B1<-bstar[[2]][[1]]
B2<-bstar[[3]][[1]]
Yhat<-B0+B1*X1+B2*X2
u<-Y-Yhat
myu<-function(u,d=1) {
n<-length(u)
c(rep(NA,d),u)[1:n]
}
u1<-myu(u)
modelu<-lm(u~u1-1)
Rho<-coef(modelu)
if(abs(Rho)>1)
break
diff<-abs(Rho-rho)
if(diff<0.00001)
break
else
rho<-Rho
}
return(coef(Model))
}
result <- boot::tsboot(A, statistic=fun, R = 1000, l=nrow(a),sim = "geom",
orig.t = TRUE)
result
Thank you.
Regards,
Lim Hock Ann
________________________________
From: Rui Barradas <ruipbarra...@sapo.pt>
To: Hock Ann Lim <lim...@yahoo.com>
Cc: R-Help <r-help@r-project.org>
Sent: Tuesday, September 18, 2012 12:03 AM
Subject: Re: [R] Problem with Stationary Bootstrap
Hello,
The problem is your function calling itself until there's no more
stack memory left.
Besides, your function doesn't make any sense. You pass an argument
'fit' then do not used it but change it's value then return itself.
Corrected:
set.seed(1)
X <- runif(10, 0, 10)
Y <- 2 + 3*X
a <- data.frame(X = X, Y = Y)
fun <- function(a){
fit <- lm(Y ~ X, data=a)
return(coef(fit))
}
result <- boot::tsboot(a, statistic = fun, R = 10, sim = "geom",
l = 10, orig.t = TRUE)
Hope this helps,
Rui Barradas
Em 17-09-2012 14:42, Hock Ann Lim escreveu:
Dear R experts,
I'm running the following stationary bootstrap programming to find the parameters estimate of a linear model:
X<-runif(10,0,10)
Y<-2+3*X
a<-data.frame(X,Y)
coef<-function(fit){
fit <- lm(Y~X,data=a)
return(coef(fit))
}
result<- tsboot(a,statistic=coef(fit),R = 10,n.sim = NROW(a),sim =
"geom",orig.t = TRUE)
Unfortunately, I got this error message from R:
Error: evaluation nested too deeply: infinite recursion / options(expressions=)?
Can someone tells me what's wrong in the programming.
Thank you.
Regards,
Lim Hock Ann [[alternative HTML version deleted]]
______________________________________________ R-help@r-project.org mailing
list 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.
______________________________________________
R-help@r-project.org mailing list
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.
