Created
April 29, 2014 22:24
-
-
Save drewlanenga/2823e9db3619bd32eac7 to your computer and use it in GitHub Desktop.
Regression technique with correlated errors, based on Cochrane-Orcutt (http://en.wikipedia.org/wiki/Cochrane%E2%80%93Orcutt_estimation)
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| ## Regression with correlated errors | |
| tolley.shuffle <- function(model,k=1e-5){ | |
| if(class(model)!="lm") | |
| stop("Object 'model' must be of class 'lm'.") | |
| uncorrelate <- function(model) { | |
| r <- residuals(model) | |
| n <- length(r) | |
| y <- model$model[,1] | |
| z <- model$model[,2:ncol(model$model)] | |
| ## Determine p and q in the ARMA(p,q) | |
| ACF <- acf(r)$acf[1:10] | |
| PACF <- pacf(r)$acf[1:10] | |
| critical <- 2/sqrt(n) | |
| a <- sum(ACF > critical) | |
| p <- sum(PACF > critical) | |
| if (a >= 4 & p <= 2) | |
| order <- c(p,0,0) | |
| else if (a <= 2 & p >= 4) | |
| order <- c(0,0,a) | |
| else order <- c(1,0,1) | |
| ## Fit the ARMA model and extract AR and MA components | |
| mod.r <- arima(r,order=order) | |
| phi <- mod.r$coef[substr(names(mod.r$coef),1,1)=='a'] | |
| theta <- mod.r$coef[substr(names(mod.r$coef),1,1)=='m'] | |
| ## Transform functions | |
| transform <- function(x,param) { | |
| if (sum(param) > 0) { | |
| x.t <- matrix(0,length(param),length(x)) | |
| for(i in 1:length(param)) { | |
| x.t[i,] <- lag(ts(x),-i)*param[i] | |
| } | |
| x.t <- colSums(x.t) | |
| } else x.t <- x | |
| return(x.t) | |
| } | |
| y.phi <- transform(y,phi) | |
| y.theta <- transform(y,theta) | |
| z.phi <- transform(z,phi) | |
| z.theta <- transform(z,theta) | |
| if (all(y.phi==y)) | |
| new.y <- y.theta | |
| else if (all(y.theta==y)) | |
| new.y <- y.phi | |
| else new.y <- y.phi/y.theta | |
| if (all(z.phi==z)) | |
| new.z <- z.theta | |
| else if (all(z.theta==z)) | |
| new.z <- z.phi | |
| else new.z <- z.phi/z.theta | |
| new.mod <- lm(new.y~new.z) | |
| return(new.mod) | |
| } | |
| ## Sets everything up for iterations until convergence | |
| mod <- model | |
| par.sums <- NULL | |
| difference <- j <- 1 | |
| while(difference > k) { | |
| j <- j+1 | |
| mod <- uncorrelate(mod) | |
| par.sums <- c(par.sums,sum(mod$coef)) | |
| difference <- par.sums[j]-par.sums[j-1] | |
| } | |
| return(list(mod,j)) | |
| } |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment