[R] Problem with Stationary Bootstrap
Rui Barradas
ruipbarradas at sapo.pt
Thu Nov 29 13:26:08 CET 2012
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 <ruipbarradas at sapo.pt>
> To: Hock Ann Lim <lim_ha at yahoo.com>
> Cc: R-Help <r-help at 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 at 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.
More information about the R-help
mailing list