[Returnanalytics-commits] r2988 - in pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code: . Covariance Matrix Integrated Regression Function Data R Tests

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Wed Sep 4 22:47:13 CEST 2013


Author: shubhanm
Date: 2013-09-04 22:47:13 +0200 (Wed, 04 Sep 2013)
New Revision: 2988

Added:
   pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/Covariance Matrix Integrated Regression Function/
   pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/Covariance Matrix Integrated Regression Function/glmi.R
   pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/Covariance Matrix Integrated Regression Function/lmi.R
   pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/Covariance Matrix Integrated Regression Function/nlsi.R
   pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/Data/
   pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/Data/Investment.csv
   pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/Data/PublicSchools.csv
   pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/Data/RealInt.csv
   pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/Data/ps.csv
   pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/R Tests/
   pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/R Tests/Cross Sectional Data.R
   pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/R Tests/HAC Data.R
   pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/R Tests/Tests.R
   pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/R Tests/Time Series Data.R
Removed:
   pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/Tests/
   pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/glmi.R
   pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/lmi.R
   pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/nlsi.R
Log:
HAC Error Test Matlab code /R code / Data 

Added: pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/Covariance Matrix Integrated Regression Function/glmi.R
===================================================================
--- pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/Covariance Matrix Integrated Regression Function/glmi.R	                        (rev 0)
+++ pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/Covariance Matrix Integrated Regression Function/glmi.R	2013-09-04 20:47:13 UTC (rev 2988)
@@ -0,0 +1,92 @@
+glmi <- function (formula, family = gaussian, data,vcov = NULL, weights, subset, 
+          na.action, start = NULL, etastart, mustart, offset, control = list(...), 
+          model = TRUE, method = "glm.fit", x = FALSE, y = TRUE, contrasts = NULL, 
+          ...) 
+{
+  call <- match.call()
+  if (is.character(family)) 
+    family <- get(family, mode = "function", envir = parent.frame())
+  if (is.function(family)) 
+    family <- family()
+  if (is.null(family$family)) {
+    print(family)
+    stop("'family' not recognized")
+  }
+  if (missing(data)) 
+    data <- environment(formula)
+  mf <- match.call(expand.dots = FALSE)
+  m <- match(c("formula", "data", "subset", "weights", "na.action", 
+               "etastart", "mustart", "offset"), names(mf), 0L)
+  mf <- mf[c(1L, m)]
+  mf$drop.unused.levels <- TRUE
+  mf[[1L]] <- as.name("model.frame")
+  mf <- eval(mf, parent.frame())
+  if (identical(method, "model.frame")) 
+    return(mf)
+  if (!is.character(method) && !is.function(method)) 
+    stop("invalid 'method' argument")
+  if (identical(method, "glm.fit")) 
+    control <- do.call("glm.control", control)
+  mt <- attr(mf, "terms")
+  Y <- model.response(mf, "any")
+  if (length(dim(Y)) == 1L) {
+    nm <- rownames(Y)
+    dim(Y) <- NULL
+    if (!is.null(nm)) 
+      names(Y) <- nm
+  }
+  X <- if (!is.empty.model(mt)) 
+    model.matrix(mt, mf, contrasts)
+  else matrix(, NROW(Y), 0L)
+  weights <- as.vector(model.weights(mf))
+  if (!is.null(weights) && !is.numeric(weights)) 
+    stop("'weights' must be a numeric vector")
+  if (!is.null(weights) && any(weights < 0)) 
+    stop("negative weights not allowed")
+  offset <- as.vector(model.offset(mf))
+  if (!is.null(offset)) {
+    if (length(offset) != NROW(Y)) 
+      stop(gettextf("number of offsets is %d should equal %d (number of observations)", 
+                    length(offset), NROW(Y)), domain = NA)
+  }
+  mustart <- model.extract(mf, "mustart")
+  etastart <- model.extract(mf, "etastart")
+  fit <- eval(call(if (is.function(method)) "method" else method, 
+                   x = X, y = Y, weights = weights, start = start, etastart = etastart, 
+                   mustart = mustart, offset = offset, family = family, 
+                   control = control, intercept = attr(mt, "intercept") > 
+                     0L))
+  if (length(offset) && attr(mt, "intercept") > 0L) {
+    fit2 <- eval(call(if (is.function(method)) "method" else method, 
+                      x = X[, "(Intercept)", drop = FALSE], y = Y, weights = weights, 
+                      offset = offset, family = family, control = control, 
+                      intercept = TRUE))
+    if (!fit2$converged) 
+      warning("fitting to calculate the null deviance did not converge -- increase 'maxit'?")
+    fit$null.deviance <- fit2$deviance
+  }
+  if (model) 
+    fit$model <- mf
+  fit$na.action <- attr(mf, "na.action")
+  if (x) 
+    fit$x <- X
+  if (!y) 
+    fit$y <- NULL
+  fit <- c(fit, list(call = call, formula = formula, terms = mt, 
+                     data = data, offset = offset, control = control, method = method, 
+                     contrasts = attr(X, "contrasts"), xlevels = .getXlevels(mt, 
+                                                                             mf)))
+  class(fit) <- c(fit$class, c("glm", "lm"))
+  fit
+  if(is.null(vcov)) {
+    se <- vcov(fit)
+  } else {
+    if (is.function(vcov))
+      se <- vcov(fit)
+    else
+      se <- vcov
+  }
+  fit = list(fit,vHaC = se) 
+  fit
+  
+}
\ No newline at end of file

Added: pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/Covariance Matrix Integrated Regression Function/lmi.R
===================================================================
--- pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/Covariance Matrix Integrated Regression Function/lmi.R	                        (rev 0)
+++ pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/Covariance Matrix Integrated Regression Function/lmi.R	2013-09-04 20:47:13 UTC (rev 2988)
@@ -0,0 +1,76 @@
+lmi <- function (formula, data,vcov = NULL, subset, weights, na.action, method = "qr", 
+          model = TRUE, x = FALSE, y = FALSE, qr = TRUE, singular.ok = TRUE, 
+          contrasts = NULL, offset, ...) 
+{
+  ret.x <- x
+  ret.y <- y
+  cl <- match.call()
+  mf <- match.call(expand.dots = FALSE)
+  m <- match(c("formula", "data", "subset", "weights", "na.action", 
+               "offset"), names(mf), 0L)
+  mf <- mf[c(1L, m)]
+  mf$drop.unused.levels <- TRUE
+  mf[[1L]] <- as.name("model.frame")
+  mf <- eval(mf, parent.frame())
+  if (method == "model.frame") 
+    return(mf)
+  else if (method != "qr") 
+    warning(gettextf("method = '%s' is not supported. Using 'qr'", 
+                     method), domain = NA)
+  mt <- attr(mf, "terms")
+  y <- model.response(mf, "numeric")
+  w <- as.vector(model.weights(mf))
+  if (!is.null(w) && !is.numeric(w)) 
+    stop("'weights' must be a numeric vector")
+  offset <- as.vector(model.offset(mf))
+  if (!is.null(offset)) {
+    if (length(offset) != NROW(y)) 
+      stop(gettextf("number of offsets is %d, should equal %d (number of observations)", 
+                    length(offset), NROW(y)), domain = NA)
+  }
+  if (is.empty.model(mt)) {
+    x <- NULL
+    z <- list(coefficients = if (is.matrix(y)) matrix(, 0, 
+                                                      3) else numeric(), residuals = y, fitted.values = 0 * 
+                y, weights = w, rank = 0L, df.residual = if (!is.null(w)) sum(w != 
+                                                                                0) else if (is.matrix(y)) nrow(y) else length(y))
+    if (!is.null(offset)) {
+      z$fitted.values <- offset
+      z$residuals <- y - offset
+    }
+  }
+  else {
+    x <- model.matrix(mt, mf, contrasts)
+    z <- if (is.null(w)) 
+      lm.fit(x, y, offset = offset, singular.ok = singular.ok, 
+             ...)
+    else lm.wfit(x, y, w, offset = offset, singular.ok = singular.ok, 
+                 ...)
+  }
+  class(z) <- c(if (is.matrix(y)) "mlm", "lm")
+  z$na.action <- attr(mf, "na.action")
+  z$offset <- offset
+  z$contrasts <- attr(x, "contrasts")
+  z$xlevels <- .getXlevels(mt, mf)
+  z$call <- cl
+  z$terms <- mt
+  if (model) 
+    z$model <- mf
+  if (ret.x) 
+    z$x <- x
+  if (ret.y) 
+    z$y <- y
+  if (!qr) 
+    z$qr <- NULL
+  #z
+  if(is.null(vcov)) {
+    se <- vcov(z)
+  } else {
+    if (is.function(vcov))
+      se <- vcov(z)
+    else
+      se <- vcov
+  }
+  z = list(z,vHaC = se) 
+  z
+}

Added: pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/Covariance Matrix Integrated Regression Function/nlsi.R
===================================================================
--- pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/Covariance Matrix Integrated Regression Function/nlsi.R	                        (rev 0)
+++ pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/Covariance Matrix Integrated Regression Function/nlsi.R	2013-09-04 20:47:13 UTC (rev 2988)
@@ -0,0 +1,178 @@
+nlsi <- function (formula, data = parent.frame(),vcov = NULL, start, control = nls.control(), 
+          algorithm = c("default", "plinear", "port"), trace = FALSE, 
+          subset, weights, na.action, model = FALSE, lower = -Inf, 
+          upper = Inf, ...) 
+{
+  formula <- as.formula(formula)
+  algorithm <- match.arg(algorithm)
+  if (!is.list(data) && !is.environment(data)) 
+    stop("'data' must be a list or an environment")
+  mf <- match.call()
+  varNames <- all.vars(formula)
+  if (length(formula) == 2L) {
+    formula[[3L]] <- formula[[2L]]
+    formula[[2L]] <- 0
+  }
+  form2 <- formula
+  form2[[2L]] <- 0
+  varNamesRHS <- all.vars(form2)
+  mWeights <- missing(weights)
+  pnames <- if (missing(start)) {
+    if (!is.null(attr(data, "parameters"))) {
+      names(attr(data, "parameters"))
+    }
+    else {
+      cll <- formula[[length(formula)]]
+      func <- get(as.character(cll[[1L]]))
+      if (!is.null(pn <- attr(func, "pnames"))) 
+        as.character(as.list(match.call(func, call = cll))[-1L][pn])
+    }
+  }
+  else names(start)
+  env <- environment(formula)
+  if (is.null(env)) 
+    env <- parent.frame()
+  if (length(pnames)) 
+    varNames <- varNames[is.na(match(varNames, pnames))]
+  lenVar <- function(var) tryCatch(length(eval(as.name(var), 
+                                               data, env)), error = function(e) -1)
+  if (length(varNames)) {
+    n <- sapply(varNames, lenVar)
+    if (any(not.there <- n == -1)) {
+      nnn <- names(n[not.there])
+      if (missing(start)) {
+        if (algorithm == "plinear") 
+          stop("no starting values specified")
+        warning("No starting values specified for some parameters.\n", 
+                "Initializing ", paste(sQuote(nnn), collapse = ", "), 
+                " to '1.'.\n", "Consider specifying 'start' or using a selfStart model", 
+                domain = NA)
+        start <- setNames(as.list(rep(1, length(nnn))), 
+                          nnn)
+        varNames <- varNames[i <- is.na(match(varNames, 
+                                              nnn))]
+        n <- n[i]
+      }
+      else stop(gettextf("parameters without starting value in 'data': %s", 
+                         paste(nnn, collapse = ", ")), domain = NA)
+    }
+  }
+  else {
+    if (length(pnames) && any((np <- sapply(pnames, lenVar)) == 
+                                -1)) {
+      message(sprintf(ngettext(sum(np == -1), "fitting parameter %s without any variables", 
+                               "fitting parameters %s without any variables"), 
+                      paste(sQuote(pnames[np == -1]), collapse = ", ")), 
+              domain = NA)
+      n <- integer()
+    }
+    else stop("no parameters to fit")
+  }
+  respLength <- length(eval(formula[[2L]], data, env))
+  if (length(n) > 0L) {
+    varIndex <- n%%respLength == 0
+    if (is.list(data) && diff(range(n[names(n) %in% names(data)])) > 
+          0) {
+      mf <- data
+      if (!missing(subset)) 
+        warning("argument 'subset' will be ignored")
+      if (!missing(na.action)) 
+        warning("argument 'na.action' will be ignored")
+      if (missing(start)) 
+        start <- getInitial(formula, mf)
+      startEnv <- new.env(hash = FALSE, parent = environment(formula))
+      for (i in names(start)) assign(i, start[[i]], envir = startEnv)
+      rhs <- eval(formula[[3L]], data, startEnv)
+      n <- NROW(rhs)
+      wts <- if (mWeights) 
+        rep(1, n)
+      else eval(substitute(weights), data, environment(formula))
+    }
+    else {
+      mf$formula <- as.formula(paste("~", paste(varNames[varIndex], 
+                                                collapse = "+")), env = environment(formula))
+      mf$start <- mf$control <- mf$algorithm <- mf$trace <- mf$model <- NULL
+      mf$lower <- mf$upper <- NULL
+      mf[[1L]] <- as.name("model.frame")
+      mf <- eval.parent(mf)
+      n <- nrow(mf)
+      mf <- as.list(mf)
+      wts <- if (!mWeights) 
+        model.weights(mf)
+      else rep(1, n)
+    }
+    if (any(wts < 0 | is.na(wts))) 
+      stop("missing or negative weights not allowed")
+  }
+  else {
+    varIndex <- logical()
+    mf <- list(0)
+    wts <- numeric()
+  }
+  if (missing(start)) 
+    start <- getInitial(formula, mf)
+  for (var in varNames[!varIndex]) mf[[var]] <- eval(as.name(var), 
+                                                     data, env)
+  varNamesRHS <- varNamesRHS[varNamesRHS %in% varNames[varIndex]]
+  m <- switch(algorithm, plinear = nlsModel.plinear(formula, 
+                                                    mf, start, wts), port = nlsModel(formula, mf, start, 
+                                                                                     wts, upper), nlsModel(formula, mf, start, wts))
+  ctrl <- nls.control()
+  if (!missing(control)) {
+    control <- as.list(control)
+    ctrl[names(control)] <- control
+  }
+  if (algorithm != "port") {
+    if (!missing(lower) || !missing(upper)) 
+      warning("upper and lower bounds ignored unless algorithm = \"port\"")
+    convInfo <- .Call(C_nls_iter, m, ctrl, trace)
+    nls.out <- list(m = m, convInfo = convInfo, data = substitute(data), 
+                    call = match.call())
+  }
+  else {
+    pfit <- nls_port_fit(m, start, lower, upper, control, 
+                         trace, give.v = TRUE)
+    iv <- pfit[["iv"]]
+    msg.nls <- port_msg(iv[1L])
+    conv <- (iv[1L] %in% 3:6)
+    if (!conv) {
+      msg <- paste("Convergence failure:", msg.nls)
+      if (ctrl$warnOnly) 
+        warning(msg)
+      else stop(msg)
+    }
+    v. <- port_get_named_v(pfit[["v"]])
+    cInfo <- list(isConv = conv, finIter = iv[31L], finTol = v.[["NREDUC"]], 
+                  nEval = c(`function` = iv[6L], gradient = iv[30L]), 
+                  stopCode = iv[1L], stopMessage = msg.nls)
+    cl <- match.call()
+    cl$lower <- lower
+    cl$upper <- upper
+    nls.out <- list(m = m, data = substitute(data), call = cl, 
+                    convInfo = cInfo, convergence = as.integer(!conv), 
+                    message = msg.nls)
+  }
+  nls.out$call$algorithm <- algorithm
+  nls.out$call$control <- ctrl
+  nls.out$call$trace <- trace
+  nls.out$na.action <- attr(mf, "na.action")
+  nls.out$dataClasses <- attr(attr(mf, "terms"), "dataClasses")[varNamesRHS]
+  if (model) 
+    nls.out$model <- mf
+  if (!mWeights) 
+    nls.out$weights <- wts
+  nls.out$control <- control
+  class(nls.out) <- "nls"
+  nls.out
+  if(is.null(vcov)) {
+    se <- vcov(nls.out)
+  } else {
+    if (is.function(vcov))
+      se <- vcov(nls.out)
+    else
+      se <- vcov
+  }
+  nls.out = list(nls.out,vHaC = se) 
+  nls.out
+  
+}
\ No newline at end of file

Added: pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/Data/Investment.csv
===================================================================
--- pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/Data/Investment.csv	                        (rev 0)
+++ pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/Data/Investment.csv	2013-09-04 20:47:13 UTC (rev 2988)
@@ -0,0 +1,21 @@
+"","GNP","Investment","Price","Interest","RealGNP","RealInv","RealInt"
+"1",596.7,90.9,0.7167,3.23,832.565927166178,126.83131017162,NA
+"2",637.7,97.4,0.7277,3.55,876.322660436993,133.84636526041,2.01518766568997
+"3",691.1,113.5,0.7436,4.04,929.397525551372,152.635825712749,1.85503366772021
+"4",756,125.7,0.7676,4.5,984.887962480459,163.757165190203,1.27245831091986
+"5",799.6,122.8,0.7906,4.19,1011.38375917025,155.325069567417,1.19364773319437
+"6",873.4,133.3,0.8254,5.16,1058.15362248607,161.497455779016,0.758279787503155
+"7",944,149.3,0.8679,5.87,1087.68291277797,172.024426777279,0.720981342379455
+"8",992.7,144.2,0.9145,5.95,1085.51120831055,157.681793329688,0.580717824634178
+"9",1077.6,166.4,0.9601,4.88,1122.3830850953,173.315279658369,-0.106331328594858
+"10",1185.9,195,1,4.5,1185.9,195,0.344182897614827
+"11",1326.4,229.8,1.0575,6.44,1254.27895981087,217.304964539007,0.68999999999999
+"12",1434.2,228.7,1.1508,7.83,1246.26346889121,198.731317344456,-0.992695035460986
+"13",1549.2,206.1,1.2579,6.25,1231.57643691867,163.844502742666,-3.05656934306569
+"14",1718,257.9,1.3234,5.5,1298.17137675684,194.87683240139,0.292908816281112
+"15",1918.3,324.1,1.4005,5.46,1369.72509817922,231.417350946091,-0.365902977180004
+"16",2163.9,386.6,1.5042,7.46,1438.57199840447,257.013694987369,0.0555016065690905
+"17",2417.8,423,1.6342,10.28,1479.50067311223,258.842246970995,1.63753224305278
+"18",2631.7,401.9,1.7842,11.77,1475.00280237642,225.255016253783,2.59119691592217
+"19",2954.1,474.9,1.9514,13.42,1513.83622014964,243.363738854156,4.0488532675709
+"20",3073,414.5,2.0688,11.02,1485.40216550657,200.357695282289,5.00380649789895

Added: pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/Data/PublicSchools.csv
===================================================================
--- pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/Data/PublicSchools.csv	                        (rev 0)
+++ pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/Data/PublicSchools.csv	2013-09-04 20:47:13 UTC (rev 2988)
@@ -0,0 +1,52 @@
+"","Expenditure","Income"
+"Alabama",275,6247
+"Alaska",821,10851
+"Arizona",339,7374
+"Arkansas",275,6183
+"California",387,8850
+"Colorado",452,8001
+"Connecticut",531,8914
+"Delaware",424,8604
+"Florida",316,7505
+"Georgia",265,6700
+"Hawaii",403,8380
+"Idaho",304,6813
+"Illinois",437,8745
+"Indiana",345,7696
+"Iowa",431,7873
+"Kansas",355,8001
+"Kentucky",260,6615
+"Louisiana",316,6640
+"Maine",327,6333
+"Maryland",427,8306
+"Massachusetts",427,8063
+"Michigan",466,8442
+"Minnesota",477,7847
+"Mississippi",259,5736
+"Missouri",274,7342
+"Montana",433,7051
+"Nebraska",294,7391
+"Nevada",359,9032
+"New Hampshire",279,7277
+"New Jersey",423,8818
+"New Mexico",388,6505
+"New York",447,8267
+"North Carolina",335,6607
+"North Dakota",311,7478
+"Ohio",322,7812
+"Oklahoma",320,6951
+"Oregon",397,7839
+"Pennsylvania",412,7733
+"Rhode Island",342,7526
+"South Carolina",315,6242
+"South Dakota",321,6841
+"Tennessee",268,6489
+"Texas",315,7697
+"Utah",417,6622
+"Vermont",353,6541
+"Virginia",356,7624
+"Washington",415,8450
+"Washington DC",428,10022
+"West Virginia",320,6456
+"Wisconsin",NA,7597
+"Wyoming",500,9096

Added: pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/Data/RealInt.csv
===================================================================
--- pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/Data/RealInt.csv	                        (rev 0)
+++ pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/Data/RealInt.csv	2013-09-04 20:47:13 UTC (rev 2988)
@@ -0,0 +1,104 @@
+"","V1"
+"1",1.99132
+"2",0.00403
+"3",2.27
+"4",0.84833
+"5",1.89112
+"6",-0.28
+"7",3.68431
+"8",1.62478
+"9",1.21599
+"10",1.28373
+"11",1.70141
+"12",3.16687
+"13",2.32779
+"14",2.26202
+"15",1.91218
+"16",3.53196
+"17",-0.31777
+"18",3.44694
+"19",1.52422
+"20",0.74268
+"21",1.31541
+"22",0.36646
+"23",3.59562
+"24",3.6574
+"25",0.97494
+"26",-0.5328
+"27",1.12681
+"28",0.31123
+"29",0.57834
+"30",1.08163
+"31",0.24977
+"32",0.23792
+"33",-0.32653
+"34",1.14733
+"35",1.19323
+"36",2.58962
+"37",0.01195
+"38",2.63842
+"39",0.42092
+"40",2.58823
+"41",-2.19809
+"42",2.89548
+"43",1.8213
+"44",0.86332
+"45",0.67496
+"46",0.34435
+"47",1.22762
+"48",-2.83991
+"49",-1.8163
+"50",-2.22965
+"51",-1.58457
+"52",-6.31184
+"53",-2.46258
+"54",-5.61479
+"55",-3.53887
+"56",1.0178
+"57",-1.58875
+"58",-1.85396
+"59",-0.2567
+"60",2.42226
+"61",-1.29502
+"62",-0.48977
+"63",1.21167
+"64",-4.85498
+"65",-3.599
+"66",0.17094
+"67",1.51603
+"68",-1.85304
+"69",-5.60478
+"70",-1.25767
+"71",0.96658
+"72",-3.09451
+"73",-5.26773
+"74",-4.03154
+"75",-1.76618
+"76",-5.73978
+"77",3.83047
+"78",0.52007
+"79",-0.18033
+"80",3.82808
+"81",3.01171
+"82",2.75289
+"83",11.74184
+"84",9.91701
+"85",3.03444
+"86",10.15143
+"87",9.29177
+"88",6.87498
+"89",2.43675
+"90",4.37202
+"91",6.77774
+"92",4.16692
+"93",5.65037
+"94",5.17734
+"95",9.41206
+"96",3.77006
+"97",4.24568
+"98",4.53154
+"99",3.39706
+"100",8.93951
+"101",4.20825
+"102",3.43461
+"103",4.30529

Added: pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/Data/ps.csv
===================================================================
--- pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/Data/ps.csv	                        (rev 0)
+++ pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/Data/ps.csv	2013-09-04 20:47:13 UTC (rev 2988)
@@ -0,0 +1,51 @@
+"","Expenditure","Income"
+"Alabama",275,0.6247
+"Alaska",821,1.0851
+"Arizona",339,0.7374
+"Arkansas",275,0.6183
+"California",387,0.885
+"Colorado",452,0.8001
+"Connecticut",531,0.8914
+"Delaware",424,0.8604
+"Florida",316,0.7505
+"Georgia",265,0.67
+"Hawaii",403,0.838
+"Idaho",304,0.6813
+"Illinois",437,0.8745
+"Indiana",345,0.7696
+"Iowa",431,0.7873
+"Kansas",355,0.8001
+"Kentucky",260,0.6615
+"Louisiana",316,0.664
+"Maine",327,0.6333
+"Maryland",427,0.8306
+"Massachusetts",427,0.8063
+"Michigan",466,0.8442
+"Minnesota",477,0.7847
+"Mississippi",259,0.5736
+"Missouri",274,0.7342
+"Montana",433,0.7051
+"Nebraska",294,0.7391
+"Nevada",359,0.9032
+"New Hampshire",279,0.7277
+"New Jersey",423,0.8818
+"New Mexico",388,0.6505
+"New York",447,0.8267
+"North Carolina",335,0.6607
+"North Dakota",311,0.7478
+"Ohio",322,0.7812
+"Oklahoma",320,0.6951
+"Oregon",397,0.7839
+"Pennsylvania",412,0.7733
+"Rhode Island",342,0.7526
+"South Carolina",315,0.6242
+"South Dakota",321,0.6841
+"Tennessee",268,0.6489
+"Texas",315,0.7697
+"Utah",417,0.6622
+"Vermont",353,0.6541
+"Virginia",356,0.7624
+"Washington",415,0.845
+"Washington DC",428,1.0022
+"West Virginia",320,0.6456
+"Wyoming",500,0.9096

Added: pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/R Tests/Cross Sectional Data.R
===================================================================
--- pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/R Tests/Cross Sectional Data.R	                        (rev 0)
+++ pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/R Tests/Cross Sectional Data.R	2013-09-04 20:47:13 UTC (rev 2988)
@@ -0,0 +1,79 @@
+library("sandwich")
+library("lmtest")
+library("strucchange")
+data("PublicSchools")
+ps <- na.omit(PublicSchools)
+ps$Income <- ps$Income * 0.0001
+fm.ps <- lm(Expenditure ~ Income + I(Income^3), data = ps)
+sqrt(diag(vcov(fm.ps)))
+sqrt(diag(vcovHC(fm.ps, type = "const")))
+sqrt(diag(vcovHC(fm.ps, type = "HC0")))
+sqrt(diag(vcovHC(fm.ps, type = "HC3")))
+sqrt(diag(vcovHC(fm.ps, type = "HC4")))
+coeftest(fm.ps, df = Inf, vcov = vcovHC(fm.ps, type = "HC0"))
+coeftest(fm.ps, df = Inf, vcov = vcovHC(fm.ps, type = "HC4"))
+plot(Expenditure ~ Income, data = ps,
+     xlab = "per capita income",
+     ylab = "per capita spending on public schools")
+inc <- seq(0.5, 1.2, by = 0.001)
+lines(inc, predict(fm.ps, data.frame(Income = inc)), col = 4, lty = 2)
+fm.ps2 <- lm(Expenditure ~ Income, data = ps)
+abline(fm.ps2, col = 4)
+text(ps[2,2], ps[2,1], rownames(ps)[2], pos = 2)
+## Willam H. Greene, Econometric Analysis, 2nd Ed.
+## Chapter 14
+## load data set, p. 385, Table 14.1
+data(PublicSchools)
+
+## omit NA in Wisconsin and scale income
+ps <- na.omit(PublicSchools)
+ps$Income <- ps$Income * 0.0001
+
+## fit quadratic regression, p. 385, Table 14.2
+fmq <- lm(Expenditure ~ Income + I(Income^2), data = ps)
+summary(fmq)
+
+## compare standard and HC0 standard errors
+## p. 391, Table 14.3
+library(sandwich)
+coef(fmq)
+sqrt(diag(vcovHC(fmq, type = "const")))
+sqrt(diag(vcovHC(fmq, type = "HC0")))
+
+if(require(lmtest)) {
+  ## compare t ratio
+  coeftest(fmq, vcov = vcovHC(fmq, type = "HC0"))
+  
+  ## White test, p. 393, Example 14.5
+  wt <- lm(residuals(fmq)^2 ~ poly(Income, 4), data = ps)
+  wt.stat <- summary(wt)$r.squared * nrow(ps)
+  c(wt.stat, pchisq(wt.stat, df = 3, lower = FALSE))
+  
+  ## Bresch-Pagan test, p. 395, Example 14.7
+  bptest(fmq, studentize = FALSE)
+  bptest(fmq)
+  
+  ## Francisco Cribari-Neto, Asymptotic Inference, CSDA 45
+  ## quasi z-tests, p. 229, Table 8
+  ## with Alaska
+  coeftest(fmq, df = Inf)[3,4]
+  coeftest(fmq, df = Inf, vcov = vcovHC(fmq, type = "HC0"))[3,4]
+  coeftest(fmq, df = Inf, vcov = vcovHC(fmq, type = "HC3"))[3,4]
+  coeftest(fmq, df = Inf, vcov = vcovHC(fmq, type = "HC4"))[3,4]
+  ## without Alaska (observation 2)
+  fmq1 <- lm(Expenditure ~ Income + I(Income^2), data = ps[-2,])
+  coeftest(fmq1, df = Inf)[3,4]
+  coeftest(fmq1, df = Inf, vcov = vcovHC(fmq1, type = "HC0"))[3,4]
+  coeftest(fmq1, df = Inf, vcov = vcovHC(fmq1, type = "HC3"))[3,4]
+  coeftest(fmq1, df = Inf, vcov = vcovHC(fmq1, type = "HC4"))[3,4]
+}
+
+## visualization, p. 230, Figure 1
+plot(Expenditure ~ Income, data = ps,
+     xlab = "per capita income",
+     ylab = "per capita spending on public schools")
+inc <- seq(0.5, 1.2, by = 0.001)
+lines(inc, predict(fmq, data.frame(Income = inc)), col = 4)
+fml <- lm(Expenditure ~ Income, data = ps)
+abline(fml)
+text(ps[2,2], ps[2,1], rownames(ps)[2], pos = 2)
\ No newline at end of file

Added: pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/R Tests/HAC Data.R
===================================================================
--- pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/R Tests/HAC Data.R	                        (rev 0)
+++ pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/R Tests/HAC Data.R	2013-09-04 20:47:13 UTC (rev 2988)
@@ -0,0 +1,17 @@
+data("RealInt")
+#OLS-based CUSUM test with quadratic spectral kernel HAC estimate:
+  ocus <- gefp(RealInt ~ 1, fit = lm, vcov = kernHAC)
+plot(ocus, aggregate = FALSE)
+sctest(ocus)
+#supF test with quadratic spectral kernel HAC estimate:
+  fs <- Fstats(RealInt ~ 1, vcov = kernHAC)
+plot(fs)
+sctest(fs)
+#Breakpoint estimation and confidence intervals with quadratic spectral kernel HAC estimate:
+  bp <- breakpoints(RealInt ~ 1)
+confint(bp, vcov = kernHAC)
+plot(bp)
+#Visualization:
+  plot(RealInt, ylab = "Real interest rate")
+lines(ts(fitted(bp), start = start(RealInt), freq = 4), col = 4)
+lines(confint(bp, vcov = kernHAC))
\ No newline at end of file

Added: pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/R Tests/Tests.R
===================================================================
--- pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/R Tests/Tests.R	                        (rev 0)
+++ pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/R Tests/Tests.R	2013-09-04 20:47:13 UTC (rev 2988)
@@ -0,0 +1,21 @@
+fpe <- read.table("http://data.princeton.edu/wws509/datasets/effort.dat")
+attach(fpe)
+lmfit = lm( change ~ setting + effort )
+sandwich(lmfit)
+Fr <- c(68,42,42,30, 37,52,24,43,
+        66,50,33,23, 47,55,23,47,
+        63,53,29,27, 57,49,19,29)
+
+Temp <- gl(2, 2, 24, labels = c("Low", "High"))
+Soft <- gl(3, 8, 24, labels = c("Hard","Medium","Soft"))
+M.user <- gl(2, 4, 24, labels = c("N", "Y"))
+Brand <- gl(2, 1, 24, labels = c("X", "M"))
+
+detg <- data.frame(Fr,Temp, Soft,M.user, Brand)
+detg.m0 <- glm(Fr ~ M.user*Temp*Soft + Brand, family = poisson, data = detg)
+summary(detg.m0)
+
+detg.mod <- glm(terms(Fr ~ M.user*Temp*Soft + Brand*M.user*Temp,
+                      keep.order = TRUE),
+                family = poisson, data = detg)
+sandwich(detg.mod)
\ No newline at end of file

Added: pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/R Tests/Time Series Data.R
===================================================================
--- pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/R Tests/Time Series Data.R	                        (rev 0)
+++ pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/R Tests/Time Series Data.R	2013-09-04 20:47:13 UTC (rev 2988)
@@ -0,0 +1,78 @@
+## Willam H. Greene, Econometric Analysis, 2nd Ed.
+## Chapter 15
+## load data set, p. 411, Table 15.1
+data(Investment)
+
+## fit linear model, p. 412, Table 15.2
+fm <- lm(RealInv ~ RealGNP + RealInt, data = Investment)
+summary(fm)
+
+## visualize residuals, p. 412, Figure 15.1
+plot(ts(residuals(fm), start = 1964),
+     type = "b", pch = 19, ylim = c(-35, 35), ylab = "Residuals")
+sigma <- sqrt(sum(residuals(fm)^2)/fm$df.residual) ## maybe used df = 26 instead of 16 ??
+abline(h = c(-2, 0, 2) * sigma, lty = 2)
+
+if(require(lmtest)) {
+  ## Newey-West covariances, Example 15.3
+  coeftest(fm, vcov = NeweyWest(fm, lag = 4))
+  ## Note, that the following is equivalent:
+  coeftest(fm, vcov = kernHAC(fm, kernel = "Bartlett", bw = 5, prewhite = FALSE, adjust = FALSE))
+  
+  ## Durbin-Watson test, p. 424, Example 15.4
+  dwtest(fm)
+  
+  ## Breusch-Godfrey test, p. 427, Example 15.6
+  bgtest(fm, order = 4)
+}
+
+## visualize fitted series
+plot(Investment[, "RealInv"], type = "b", pch = 19, ylab = "Real investment")
+lines(ts(fitted(fm), start = 1964), col = 4)
+
+## 3-d visualization of fitted model
+if(require(scatterplot3d)) {
+  s3d <- scatterplot3d(Investment[,c(5,7,6)],
+                       type = "b", angle = 65, scale.y = 1, pch = 16)
+  s3d$plane3d(fm, lty.box = "solid", col = 4)
+}
+## fit investment equation
+data(Investment)
+fm <- lm(RealInv ~ RealGNP + RealInt, data = Investment)
+
+## Newey & West (1994) compute this type of estimator
+NeweyWest(fm)
+
+## The Newey & West (1987) estimator requires specification
+## of the lag and suppression of prewhitening
+NeweyWest(fm, lag = 4, prewhite = FALSE)
+
+## bwNeweyWest() can also be passed to kernHAC(), e.g.
+## for the quadratic spectral kernel
+kernHAC(fm, bw = bwNeweyWest)
+
+curve(kweights(x, kernel = "Quadratic", normalize = TRUE),
+      from = 0, to = 3.2, xlab = "x", ylab = "k(x)")
+curve(kweights(x, kernel = "Bartlett", normalize = TRUE),
+      from = 0, to = 3.2, col = 2, add = TRUE)
+curve(kweights(x, kernel = "Parzen", normalize = TRUE),
+      from = 0, to = 3.2, col = 3, add = TRUE)
+curve(kweights(x, kernel = "Tukey", normalize = TRUE),
+      from = 0, to = 3.2, col = 4, add = TRUE)
+curve(kweights(x, kernel = "Truncated", normalize = TRUE),
+      from = 0, to = 3.2, col = 5, add = TRUE)
+
+## fit investment equation
+data(Investment)
+fm <- lm(RealInv ~ RealGNP + RealInt, data = Investment)
+
+## compute quadratic spectral kernel HAC estimator
+kernHAC(fm)
+kernHAC(fm, verbose = TRUE)
+
+## use Parzen kernel instead, VAR(2) prewhitening, no finite sample
+## adjustment and Newey & West (1994) bandwidth selection
+kernHAC(fm, kernel = "Parzen", prewhite = 2, adjust = FALSE,
+        bw = bwNeweyWest, verbose = TRUE)
+## compare with estimate under assumption of spheric errors
+vcov(fm)
\ No newline at end of file

Deleted: pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/glmi.R
===================================================================
--- pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/glmi.R	2013-09-04 19:09:13 UTC (rev 2987)
+++ pkg/PerformanceAnalytics/sandbox/Shubhankit/Week6-7/Code/glmi.R	2013-09-04 20:47:13 UTC (rev 2988)
@@ -1,92 +0,0 @@
-glmi <- function (formula, family = gaussian, data,vcov = NULL, weights, subset, 
-          na.action, start = NULL, etastart, mustart, offset, control = list(...), 
-          model = TRUE, method = "glm.fit", x = FALSE, y = TRUE, contrasts = NULL, 
-          ...) 
-{
-  call <- match.call()
-  if (is.character(family)) 
-    family <- get(family, mode = "function", envir = parent.frame())
-  if (is.function(family)) 
-    family <- family()
-  if (is.null(family$family)) {
-    print(family)
-    stop("'family' not recognized")
-  }
-  if (missing(data)) 
-    data <- environment(formula)
-  mf <- match.call(expand.dots = FALSE)
-  m <- match(c("formula", "data", "subset", "weights", "na.action", 
-               "etastart", "mustart", "offset"), names(mf), 0L)
-  mf <- mf[c(1L, m)]
-  mf$drop.unused.levels <- TRUE
-  mf[[1L]] <- as.name("model.frame")
-  mf <- eval(mf, parent.frame())
-  if (identical(method, "model.frame")) 
-    return(mf)
-  if (!is.character(method) && !is.function(method)) 
-    stop("invalid 'method' argument")
-  if (identical(method, "glm.fit")) 
-    control <- do.call("glm.control", control)
-  mt <- attr(mf, "terms")
-  Y <- model.response(mf, "any")
-  if (length(dim(Y)) == 1L) {
-    nm <- rownames(Y)
-    dim(Y) <- NULL
-    if (!is.null(nm)) 
-      names(Y) <- nm
-  }
-  X <- if (!is.empty.model(mt)) 
-    model.matrix(mt, mf, contrasts)
-  else matrix(, NROW(Y), 0L)
-  weights <- as.vector(model.weights(mf))
-  if (!is.null(weights) && !is.numeric(weights)) 
-    stop("'weights' must be a numeric vector")
-  if (!is.null(weights) && any(weights < 0)) 
-    stop("negative weights not allowed")
-  offset <- as.vector(model.offset(mf))
-  if (!is.null(offset)) {
-    if (length(offset) != NROW(Y)) 
-      stop(gettextf("number of offsets is %d should equal %d (number of observations)", 
-                    length(offset), NROW(Y)), domain = NA)
-  }
-  mustart <- model.extract(mf, "mustart")
-  etastart <- model.extract(mf, "etastart")
-  fit <- eval(call(if (is.function(method)) "method" else method, 
-                   x = X, y = Y, weights = weights, start = start, etastart = etastart, 
-                   mustart = mustart, offset = offset, family = family, 
-                   control = control, intercept = attr(mt, "intercept") > 
-                     0L))
-  if (length(offset) && attr(mt, "intercept") > 0L) {
-    fit2 <- eval(call(if (is.function(method)) "method" else method, 
-                      x = X[, "(Intercept)", drop = FALSE], y = Y, weights = weights, 
-                      offset = offset, family = family, control = control, 
-                      intercept = TRUE))
-    if (!fit2$converged) 
-      warning("fitting to calculate the null deviance did not converge -- increase 'maxit'?")
-    fit$null.deviance <- fit2$deviance
-  }
-  if (model) 
-    fit$model <- mf
-  fit$na.action <- attr(mf, "na.action")
-  if (x) 
-    fit$x <- X
-  if (!y) 
-    fit$y <- NULL
-  fit <- c(fit, list(call = call, formula = formula, terms = mt, 
[TRUNCATED]

To get the complete diff run:
    svnlook diff /svnroot/returnanalytics -r 2988


More information about the Returnanalytics-commits mailing list