# R/ols.R In gets: General-to-Specific (GETS) Modelling and Indicator Saturation Methods

#### Documented in ols

```ols <-
function(y, x, untransformed.residuals=NULL, tol=1e-07,
LAPACK=FALSE, method=3, user.fun=NULL, user.options=NULL)
{

##to do number 1:
## - in version 0.14: remove the method=0 option
## - rename ols to estFun?
## - merge user.fun and user.options into a single argument,
## user.estimator, which is a list containing at least one
## entry, name, i.e. the name of the estimator-function
##
##to do number 2:
## - split estFun into two functions, estFun and vcovFun

##user-specified:
if(method==0){
message("NOTE: method=0 will be deprecated in future versions")
user.options <- c(list(y = y, x = x), user.options)
out <- do.call(user.fun, user.options, envir = .GlobalEnv)
}

##fastest (usually only for estimates):
if(method==1){
out <- list()
qx <- qr(x, tol, LAPACK=LAPACK)
out <- c(out, qx)
out\$coefficients <- solve.qr(qx, y, tol=tol)
}

##second fastest (slightly more output):
if(method==2){
out <- list()
qx <- qr(x, tol, LAPACK=LAPACK) ## compute qr-decomposition of x
out <- c(out, qx)
out\$coefficients <- solve.qr(qx, y, tol=tol)
out\$xtxinv <- chol2inv(qx\$qr, LINPACK=FALSE) #(x'x)^-1
out\$fit <- as.vector(x %*% out\$coefficients)
out\$residuals <- y - out\$fit
}

##ordinary vcov:
if(method==3){
out <- list()
out\$n <- length(y)
if(is.null(x)){ out\$k <- 0 }else{ out\$k <- NCOL(x) }
out\$df <- out\$n - out\$k
if(out\$k > 0){
qx <- qr(x, tol, LAPACK=LAPACK) ## compute qr-decomposition of x
out <- c(out, qx)
out\$coefficients <- solve.qr(qx, y, tol=tol)
out\$xtxinv <- chol2inv(qx\$qr, LINPACK=FALSE) #(x'x)^-1
out\$fit <- as.vector(x %*% out\$coefficients)
}else{
out\$fit <- rep(0, out\$n)
}
out\$residuals <- y - out\$fit
out\$residuals2 <- out\$residuals^2
if(out\$k>0){
out\$vcov <- out\$sigma2 * out\$xtxinv
}
}

##white vcov:
if(method==4){
out <- list()
out\$n <- length(y)
if(is.null(x)){ out\$k <- 0 }else{ out\$k <- NCOL(x) }
out\$df <- out\$n - out\$k
if(out\$k > 0){
qx <- qr(x, tol, LAPACK=LAPACK) ## compute qr-decomposition of x
out <- c(out, qx)
out\$coefficients <- solve.qr(qx, y, tol=tol)
out\$xtxinv <- chol2inv(qx\$qr, LINPACK=FALSE) #(x'x)^-1
out\$fit <- as.vector(x %*% out\$coefficients)
}else{
out\$fit <- rep(0, out\$n)
}
out\$residuals <- y - out\$fit
out\$residuals2 <- out\$residuals^2
if(out\$k>0){
out\$omegahat <- crossprod(x, x*out\$residuals2)
out\$vcov <- out\$xtxinv %*% out\$omegahat %*% out\$xtxinv
}
}

##newey-west vcov:
if(method==5){
out <- list()
out\$n <- length(y)
if(is.null(x)){ out\$k <- 0 }else{ out\$k <- NCOL(x) }
out\$df <- out\$n - out\$k
if(out\$k>0){
qx <- qr(x, tol, LAPACK=LAPACK) ## compute qr-decomposition of x
out <- c(out, qx)
out\$coefficients <- solve.qr(qx, y, tol=tol)
out\$xtxinv <- chol2inv(qx\$qr, LINPACK=FALSE) #(x'x)^-1
out\$fit <- as.vector(x %*% out\$coefficients)
}else{
out\$fit <- rep(0, out\$n)
}
out\$residuals <- y - out\$fit
out\$residuals2 <- out\$residuals^2

if(out\$k>0){
iL <- round(out\$n^(1/4), digits=0)
vW <- 1 - 1:iL/(iL+1)
vWsqrt <- sqrt(vW)

mSum <- 0
for(l in 1:iL){
}

out\$omegahat <- mS0 + mSum
out\$vcov <- out\$xtxinv %*% out\$omegahat %*% out\$xtxinv
} #end if(out\$k>0)

}

##log-variance w/ordinary vcov (note: y = log(e^2)):
if(method==6){
out <- list()
out\$n <- length(y)
if(is.null(x)){ out\$k <- 0 }else{ out\$k <- NCOL(x) }
out\$df <- out\$n - out\$k
if(out\$k > 0){
qx <- qr(x, tol, LAPACK=LAPACK) ## compute qr-decomposition of x
out <- c(out, qx)
out\$coefficients <- solve.qr(qx, y, tol=tol)
out\$xtxinv <- chol2inv(qx\$qr, LINPACK=FALSE) #(x'x)^-1
out\$fit <- as.vector(x %*% out\$coefficients)
}else{
out\$fit <- rep(0, out\$n)
}
out\$residuals <- y - out\$fit #residuals of AR-X representation
out\$residuals2 <- out\$residuals^2
if(out\$k>0){
out\$vcov <- out\$sigma2 * out\$xtxinv
}
##log-variance part:
out\$Elnz2 <- -log(mean(exp(out\$residuals)))
out\$var.fit <- exp(out\$fit - out\$Elnz2)
out\$std.residuals <- untransformed.residuals/sqrt(out\$var.fit)
out\$logl <- -out\$n*log(2*pi)/2 - sum(log(out\$var.fit))/2 - sum(untransformed.residuals^2/out\$var.fit)/2
}

##result:
return(out)

}
```

## Try the gets package in your browser

Any scripts or data that you put into this service are public.

gets documentation built on April 4, 2018, 5:03 p.m.