[Returnanalytics-commits] r3295 - in pkg/FactorAnalytics: R sandbox

Mon Dec 23 07:14:36 CET 2013

Author: chenyian
Date: 2013-12-23 07:14:36 +0100 (Mon, 23 Dec 2013)
New Revision: 3295

Modified:
   pkg/FactorAnalytics/R/.Rhistory
   pkg/FactorAnalytics/R/fitFundamentalFactorModel.R
   pkg/FactorAnalytics/sandbox/test.vignette.r
Log:
debug: 1. assign the right colnames to industry model in fitFundamentalFactorModel.R
2. assign the correct weights to wls method in fitFundamentalFactorModel.R

Modified: pkg/FactorAnalytics/R/.Rhistory
===================================================================

--- pkg/FactorAnalytics/R/.Rhistory	2013-12-20 01:54:08 UTC (rev 3294)
+++ pkg/FactorAnalytics/R/.Rhistory	2013-12-23 06:14:36 UTC (rev 3295)
@@ -1,512 +1,512 @@
-eigenvector <- eigen(cov(data))$vectors
-eigenvalues <- eigen(cov(data))$values
-abline(a=0,b=eigenvector[2,1]/eigenvector[1,1],col="red")
-abline(a=0,b=eigenvector[2,2]/eigenvector[1,2],col="red")
-covvar <- cbind(c(2,1),c(1,1))
-data <- rmvnorm(100,mean=c(0,0),sigma=covvar)
-eigenvector <- eigen(cov(data))$vectors
-eigenvector
-eigenvalues <- eigen(cov(data))$values
-eigenvalues
-plot(data)
-abline(a=0,b=eigenvector[2,1]/eigenvector[1,1],col="red")
-abline(a=0,b=eigenvector[2,2]/eigenvector[1,2],col="red")
-eigenvector <- eigen(cor(data))$vectors
-eigenvector
-eigenvalues <- eigen(cor(data))$values
-eigenvalues
-plot(data)
-abline(a=0,b=eigenvector[2,1]/eigenvector[1,1],col="red")
-abline(a=0,b=eigenvector[2,2]/eigenvector[1,2],col="red")
-covvar <- cbind(c(2,-1),c(-1,1))
-data <- rmvnorm(100,mean=c(0,0),sigma=covvar)
-eigenvector <- eigen(cor(data))$vectors
-eigenvector
-eigenvalues <- eigen(cor(data))$values
-eigenvalues
-plot(data)
-abline(a=0,b=eigenvector[2,1]/eigenvector[1,1],col="red")
-abline(a=0,b=eigenvector[2,2]/eigenvector[1,2],col="red")
-covvar <- cbind(c(0,-1),c(-1,0))
-data <- rmvnorm(100,mean=c(0,0),sigma=covvar)
-covvar <- cbind(c(0,1),c(1,0))
-data <- rmvnorm(100,mean=c(0,0),sigma=covvar)
-eigen(covvar)
-covvar <- cbind(c(0,1,0),c(1,0,1),c(0,1,0))
-covvar
-eigen(covvar)
-covvar <- cbind(c(1,1,0),c(1,1,1),c(0,1,1))
-covvar
-eigen(covvar)
-covvar <- cbind(c(2,1,0),c(1,1,1),c(0,1,1))
-covvar
-eigen(covvar)
-cor(covvar)
-covvar <- cbind(c(1,1,0),c(1,1,1),c(0,1,1))
-covvar
-eigen(covvar)
-cor(covvar)
-chol <- cbind(c(1,p,q),c(0,1,r),c(0,0,1))
-covvar <- chol%*%t(chol)
-covvar
-p <- 0.9
-q <- 0.1
-r <- 0.8
-chol <- cbind(c(1,p,q),c(0,1,r),c(0,0,1))
-covvar <- chol%*%t(chol)
-covvar
-is.positive.definite(chol)
-eigen(covvar)
-covvar
-chol <- cbind(c(2,p,q),c(0,1,r),c(0,0,1))
-covvar <- chol%*%t(chol)
-covvar
-is.positive.definite(chol)
-eigen(covvar)
-chol <- cbind(c(20,p,q),c(0,1,r),c(0,0,1))
-covvar <- chol%*%t(chol)
-covvar
-eigen(covvar)
-data <- rmvnorm(100,mean=c(0,0,0),sigma=covvar)
-install.packages("scatterplot3d")
-library(mvtnorm)
-library(scatterplot3d)
-? scatterplot3d
-scatterplot3d(data)
-trans3d(data)
-trans3d(data[,1],data[,2],data[,3])
-scatterplot3d(data,highlight.3d=TRUE, col.axis="blue",
-col.grid="lightblue", main="scatterplot3d - 1", pch=20)
-eigen(covvar)
-eigenvector <- eigen(cor(data))$vectors
-eigenvector
-eigenvector <- eigen(cov(data))$vectors
-eigenvector
-p <- 10
-chol <- cbind(c(1,p,q),c(0,1,r),c(0,0,1))
-covvar <- chol%*%t(chol)
-covvar
-p <- 1
-chol <- cbind(c(1,p,q),c(0,1,r),c(0,0,1))
-covvar <- chol%*%t(chol)
-covvar
-chol <- cbind(c(1,p,q),c(0,0,r),c(0,0,1))
-covvar <- chol%*%t(chol)
-covvar
-p <- .1
-chol <- cbind(c(1,p,q),c(0,0,r),c(0,0,1))
-covvar <- chol%*%t(chol)
-p <- .1
-p <- .1
-q <- 0.1
-r <- 0.8
-chol <- cbind(c(1,p,q),c(0,0,r),c(0,0,1))
-covvar <- chol%*%t(chol)
-covvar
-chol <- cbind(c(1,p,q),c(0,1,r),c(0,0,1))
-covvar <- chol%*%t(chol)
-covvar
-chol <- cbind(c(1,p,q),c(0,1,r),c(0,0,0.7))
-covvar <- chol%*%t(chol)
-covvar
-chol <- cbind(c(1,p,q),c(0,1,r),c(0,0,.65))
-covvar <- chol%*%t(chol)
-covvar
-is.positive.definite(chol)
-eigen(covvar)
-r <- 10
-chol <- cbind(c(1,p,q),c(0,1,r),c(0,0,.65))
-covvar <- chol%*%t(chol)
-covvar
-r <- 1
-chol <- cbind(c(1,p,q),c(0,1,r),c(0,0,.65))
-covvar <- chol%*%t(chol)
-covvar
-r <- 0.1
-chol <- cbind(c(1,p,q),c(0,1,r),c(0,0,.65))
-covvar <- chol%*%t(chol)
-covvar
-chol <- cbind(c(1,p,q),c(0,1,r),c(0,0,.9))
-covvar <- chol%*%t(chol)
-covvar
-chol <- cbind(c(1,p,q),c(0,1,r),c(0,0,.98))
-covvar <- chol%*%t(chol)
-covvar
-chol <- cbind(c(1,p,q),c(0,1,r),c(0,0,.999))
-covvar <- chol%*%t(chol)
-covvar
-eigen(covvar)
-chol <- cbind(c(1,p,q),c(0,1,r),c(0,0,.997))
-covvar <- chol%*%t(chol)
-covvar
-chol <- cbind(c(1,p,q),c(0,1,r),c(0,0,.9999))
-covvar <- chol%*%t(chol)
-covvar
-r <- 0.5
-chol <- cbind(c(1,p,q),c(0,1,r),c(0,0,.9999))
-covvar <- chol%*%t(chol)
-covvar
-chol <- cbind(c(1,p,q),c(0,1,r),c(0,0,.8))
-covvar <- chol%*%t(chol)
-covvar
-chol <- cbind(c(1,p,q),c(0,1,r),c(0,0,.9))
-covvar <- chol%*%t(chol)
-covvar
-eigen(covvar)
-chol <- cbind(c(1,p,q),c(0,1,r),c(0,0,.9))
-data <- rmvnorm(100,mean=c(0,0,0),sigma=covvar)
-pc <- princomp(data)
-summary(pc)
-plot(pc)
-loadings(pc)
-eigen(cov(data))
-data <- rmvnorm(100,mean=c(1,0,0),sigma=covvar)
-eigen(cov(data))
-pc <- princomp(data)
-summary(pc)
-loadings(pc)
-data <- rmvnorm(100,mean=c(0,0,0),sigma=covvar)
-# download data small scale experiment
-# using the finance sector provided by CNN money
-library(quantmod)
-library(PerformanceAnalytics)
-symbol.vec=c("AAN","AB","ACAS","ACY","AFL","AIG","AMG","AXP","BAC","BGCP",
-"C","CCNE","DB","GS","HCC","IHC","JPM","KEY","PLFE","TCHC")
-getSymbols(symbol.vec, from ="2000-01-03", to = "2012-05-10")
-# extract monthly adjusted closing prices
-l <- length(symbol.vec)
-db.m.price <- to.monthly(AAN)[, "AAN.Adjusted", drop=FALSE]
-colnames(db.m.price) <- "AAN"
-db.m.ret <- CalculateReturns(db.m.price, method="compound")[-1,]
-for (i in (2:l)) {
-name.price <-  paste(symbol.vec[i],"m","price",sep=".")
-stock <- as.name(symbol.vec[i])
-db.m.new <- to.monthly(eval(stock))[,"eval(stock).Adjusted",drop=FALSE]
-colnames(db.m.new) <- symbol.vec[i]
-db.m.price <- cbind(db.m.price,db.m.new)
-# calculate log-returns
-db.m.ret.new <- CalculateReturns(db.m.new, method="compound")[-1,]
-db.m.ret <- cbind(db.m.ret,db.m.ret.new)
+nrow=6)
+gplot(t(adj.mat),gmode="digraph",label=c(1,2,3,4,5,6),vertex.cex=2,
+arrowhead.cex = 1)
+gplot(t(adj.mat1),gmode="digraph",label=c(1,2,3,4,5,6),vertex.cex=2,
+arrowhead.cex = 1)
+gplot(t(adj.mat1),gmode="digraph",label=c(1,2,3,4,5,6),vertex.cex=2,
+arrowhead.cex = 1)
+gplot(t(adj.mat),gmode="digraph",label=c(1,2,3,4,5,6),vertex.cex=2,
+arrowhead.cex = 1)
+rm(list=ls())
+library(factorAnalytics)
+library(fEcofin)
+ts.berndt<-xts(berndtInvest[,-1],as.Date(berndtInvest[,1]))
+data(stat.fm.data)
+View(sfm.dat)
+install.packages("~/R-project/returnanalytics/pkg/FactorAnalytics/sandbox/fEcofin_2100.77.zip-UWunsafe", repos = NULL)
+install.packages("~/R-project/returnanalytics/pkg/FactorAnalytics/sandbox/fEcofin_2100.77.zip", repos = NULL)
+library(fEcofin)
+? fEconfin
+fEconfin
+help(pakcage=fEconfin)
+help(package=fEconfin)
+help(package=fEcofin)
+data(berndtInvest)
+ts.berndt<-xts(berndtInvest[,-1],as.Date(berndtInvest[,1]))
+berndt<-ts.berndt['1978/1987']
+returns<-berndt[,c(-10,-17)]
+tickers <- names(returns)
+num.tickers <- length(tickers)
+dates <- index(returns)
+num.dates <- length(dates)
+sector<-c("OTHER","OTHER","OIL","TECH","TECH","OIL","OTHER","OTHER",
+"TECH","OIL","OIL","OTHER","TECH","OIL","OTHER")
+stacked.returns <- data.frame(
+DATE=rep(dates,num.tickers),
+TICKER=rep(tickers,each=num.dates),
+RETURN=c(coredata(returns)),
+SECTOR=rep(sector,each=num.dates),
+stringsAsFactors=FALSE)
+head(stacked.returns)
+barra <- fitFundamentalFactorModel(data=stacked.returns,exposure.names="SECTOR",
+datevar="DATE",returnsvar="RETURN",
+assetvar="TICKER",wls=TRUE)
+head(barra$factor.returns)
+barra$beta
+barra.cov <- factorModelCovariance(barra$beta, barra$factor.cov$cov, barra$resid.variance)
+barra.cor <- cov2cor(barra.cov)
+round(barra.cor,2)
+View(stacked.returns)
+returns = berndtInvest[,-c(1,11,18)]
+n.stocks = ncol(returns)
+tech.dum = oil.dum = other.dum = matrix(0,n.stocks,1)
+tech.dum[c(4,5,9,13),] = 1
+oil.dum[c(3,6,10,11,14),] = 1
+other.dum = 1 - tech.dum - oil.dum
+B = cbind(tech.dum,oil.dum,other.dum)
+dimnames(B) = list(colnames(returns),c("TECH","OIL","OTHER"))
+B
+returns = t(returns)
+F.hat = solve(crossprod(B))%*%t(B)%*%returns
+E.hat = returns - B%*%F.hat
+diagD.hat = apply(E.hat,1,var)
+Dinv.hat = diag(diagD.hat^(-1))
+H = solve(t(B)%*%Dinv.hat%*%B)%*%t(B)%*%Dinv.hat
+round(H[,1:8],5)
+apply(H,1,sum)
+F.hat = H%*%returns
+E.hat = returns - B%*%F.hat
+diagD.hat = apply(E.hat,1,var)
+F.hat = t(F.hat)
+F.hat
+round(H[,1:8],5)
+H
+B
+returns
+View(berndtInvest)
+H
+names(barra)
+barra$factor.returns
+H
+F.hat
+head(barra$factor.returns)
+head(F.hat)
+oil.f.ret <- barra$returns[,1]
+oth.f.ret <- barra$returns[,1]+barra$returns[,2]
+tech.f.ret <- barra$returns[,1]+barra$returns[,3]
+head(oil.f.ret)
+oil.f.ret <- barra$factor.returns[,1]
+tech.f.ret <- barra$factor.returns[,1]+ barra$factor.returns[,3]
+oth.f.ret <- barra$factor.returns[,1]+barra$factor.returns[,2]
+head(oil.f.ret)
+head(F.hat)
+head(F.hat)[,2]
+head(barra$factor.returns)
+head(F.hat)
+head(cbind(tech.f.ret,oil.f.ret,oth.f.ret))
+tickers
+cbind(tickers,sector)
+barra$beta
+names(barra)
+? fitFundamentalFactorModel
+barra.cov2 <- barra$factor.cov
+barra.cov2
+barra.cov2 <- barra$returns.cov
+barra.cov2
+barra.cor2 <- cov2cor(barra.cov2)
+barra.cor2 <- cov2cor(barra.cov2$cov)
+cov.ind = B%*%var(F.hat)%*%t(B) + diag(diagD.hat)
+sd = sqrt(diag(cov.ind))
+cor.ind = cov.ind/outer(sd,sd)
+cor.samp <- cor(t(returns))
+View(cor.samp)
+View(barra.cor2)
+data(berndtInvest)
+ts.berndt<-xts(berndtInvest[,-1],as.Date(berndtInvest[,1]))
+berndt<-ts.berndt['1978/1987']
+returns<-berndt[,c(-10,-17)]
+tickers <- names(returns)
+num.tickers <- length(tickers)
+dates <- index(returns)
+num.dates <- length(dates)
+sector<-c("OTHER","OTHER","OIL","TECH","TECH","OIL","OTHER","OTHER",
+"TECH","OIL","OIL","OTHER","TECH","OIL","OTHER")
+stacked.returns <- data.frame(
+DATE=rep(dates,num.tickers),
+TICKER=rep(tickers,each=num.dates),
+RETURN=c(coredata(returns)),
+SECTOR=rep(sector,each=num.dates),
+stringsAsFactors=FALSE)
+head(stacked.returns)
+data=stacked.returns
+exposure.names="SECTOR"
+datevar="DATE",
+datevar="DATE"
+returnsvar="RETURN"
+assetvar="TICKER"
+wls=TRUE
+full.resid.cov=FALSE
+assets = unique(data[[assetvar]])
+timedates = as.Date(unique(data[[datevar]]))
+data[[datevar]] <- as.Date(data[[datevar]])
+if (length(timedates) < 2)
+stop("At least two time points, t and t-1, are needed for fitting the factor model.")
+if (!is(exposure.names, "vector") || !is.character(exposure.names))
+stop("exposure argument invalid---must be character vector.")
+if (!is(assets, "vector") || !is.character(assets))
+stop("assets argument invalid---must be character vector.")
+wls <- as.logical(wls)
+full.resid.cov <- as.logical(full.resid.cov)
+robust.scale = FALSE
+standardized.factor.exposure = FALSE
+numTimePoints <- length(timedates)
+numExposures <- length(exposure.names)
+numAssets <- length(assets)
+# check if exposure.names are numeric, if not, create exposures. factors by dummy variables
+which.numeric <- sapply(data[, exposure.names, drop = FALSE],is.numeric)
+exposures.numeric <- exposure.names[which.numeric]
+# industry factor model
+exposures.factor <- exposure.names[!which.numeric]
+if (length(exposures.factor) > 1) {
+stop("Only one nonnumeric variable can be used at this time.")
 }
-head(db.m.price)
-dim(db.m.price)
-dim(db.m.ret)
-corr.m <- cor(db.m.ret)
-corr.m.inv <- solve(corr.m)
-db.pc <- princomp(db.m.ret)
-summary(db.pc)
-centrality <- loadings(db.pc)[,1]
-centrality
-eigen(corr.m.inv)
-centrality
-centrality.inv <- eigen(corr.m.inv)$vectors[,1]
-eigen(corr.m)$vectors[,1]
-cov.m.inv <- solve(cov(db.m.ret))
-centrality.inv <- eigen(cov.m.inv)$vectors[,1]
-centrality.inv
-eigen(cov(db.m.ret))$vectors[,1]
-centrality
-centrality.inv
-head(db.m.ret)
-names(centrality.inv) <- colnames(db.m.ret)
-centrality.inv
-covvar <- cbind(c(0.8,0.1,0.1),c(0.8,0.1,.1),c(.8,.1,.1))
-covvar
-covvar <- cbind(c(0.8,0.8,0.8),c(0.1,0.1,.1),c(.1,.1,.1))
-covvar
-eigen(covvar)
-covvar <- cbind(c(0.5,0.5,0.5),c(0.4,0.4,.4),c(.1,.1,.1))
-covvar
-eigen(covvar)
-sum(eigen(covvar)$vectors[,1]^2)
-eigen(covvar)$vectors[,1]^2
-sd(eigen(covvar)$vectors[,1])
-covvar <- cbind(replicate(3,c(.33,.33,.33))
-covvar <- cbind(replicate(3,c(.33,.33,.33)))
-covvar
-covvar <- cbind(replicate(3,c(.33,.33,.33)))
-covvar
-eigen(covvar)
-sd(eigen(covvar)$vectors[,1])
-covvar <- cbind(c(0.5,0.4,0.5),c(0.4,0.5,.4),c(.1,.1,.1))
-covvar
-eigen(covvar)
-covvar <- cbind(c(0.5,0.4,0.3),c(0.4,0.5,.6),c(.1,.1,.1))
-covvar
-eigen(covvar)
-sd(eigen(covvar)$vectors[,1])
-covvar <- cbind(c(0.5,0.4,0.3),c(0.4,0.5,.5),c(.1,.1,.2))
-covvar
-eigen(covvar)
-sd(eigen(covvar)$vectors[,1])
-covvar <- cbind(c(0,0.7,0.5),c(0.9,0,.5),c(.1,0.3,0))
-covvar
-eigen(covvar)
-eigen(t(covvar)
-eigen(t(covvar))
-eigen(t(covvar))
-covvar
-t(covvar)
-eigen(t(covvar))
-covvar <- cbind(c(0.1,0.7,0.5),c(0.8,0.1,.2),c(.1,0.2,0.3))
-covvar
-t(covvar)
-eigen(t(covvar))
-sd(eigen(covvar)$vectors[,1])
-eigen(covvar)
-t(covvar)
-eigen(t(covvar))
-covvar <- cbind(c(0.8,0.8,0.8),c(0.1,0.1,.1),c(.1,0.1,0.1))
-covvar
-eigen(covvar)
-t(covvar)
-eigen(t(covvar))
-sd(eigen(covvar)$vectors[,1])
-sd(eigen(t(covvar))$vectors[,1])
-t(covvar)
-library("rmgarch")
-? dcc
-? DCC.fit
-? dccfit
-eigen(diag(2))
-eigen(matrix(rep(1,4),nrow=2))
-eigen(matrix(c(1,-1,-1,1),nrow=2))
-covvar <- matrix(rep(1,9),nrow=3)
-n <- length(covvar[1,])
-alpha <- eigen(cov(covvar))$values[1] -10^(-3)
-solve(diag(n)-alpha*cov(covvar))%*%rep(1,n)
-alpha <- eigen(covvar)$values[1] -10^(-3)
-solve(diag(n)-alpha*cov(covvar))%*%rep(1,n)
-eigen(covvar)$values
-alpha <- eigen(covvar)$values[1] -10^(-3)
-solve(diag(n)-alpha*cov(covvar))%*%rep(1,n)
-kalz <- function(covvar) {
-n <- length(covvar[1,])
-alpha <- eigen(covvar)$values[1] -10^(-3)
-kalz.ec <- solve(diag(n)-alpha*cov(covvar))%*%rep(1,n)
-return(kalz.ec)
+exposures.factor
+regression.formula <- paste("~", paste(exposure.names, collapse = "+"))
+if (length(exposures.factor)) {
+regression.formula <- paste(regression.formula, "- 1")
+data[, exposures.factor] <- as.factor(data[,exposures.factor])
+exposuresToRecode <- names(data[, exposure.names, drop = FALSE])[!which.numeric]
+contrasts.list <- lapply(seq(length(exposuresToRecode)),
+function(i) function(n, m) contr.treatment(n, contrasts = FALSE))
+names(contrasts.list) <- exposuresToRecode
+}    else {
+contrasts.list <- NULL
 }
-kalz(covvar)
-covvar2 <- matrix(rep(.1,9),nrow=3)
-kalz(covvar2)
-kalz <- function(covvar) {
-n <- length(covvar[1,])
-alpha <- eigen(covvar)$values[1] -10^(-1)
-kalz.ec <- solve(diag(n)-alpha*cov(covvar))%*%rep(1,n)
-return(kalz.ec)
+# turn characters into formula
+regression.formula <- eval(parse(text = paste(returnsvar,regression.formula)))
+# RETURN ~ BOOK2MARKET
+regression.formula
+wls.classic <- function(xdf, modelterms, conlist, w) {
+assign("w", w, pos = 1)
+model <- try(lm(formula = modelterms, data = xdf, contrasts = conlist,
+weights = w, singular.ok = FALSE))
+if (is(model, "Error")) {
+mess <- geterrmessage()
+nn <- regexpr("computed fit is singular", mess)
+if (nn > 0) {
+cat("At time:", substring(mess, nn), "\n")
+model <- lm(formula = modelterms, data = xdf,
+contrasts = conlist, weights = w)
 }
-# example of all 1 matrix
-covvar <- matrix(rep(1,9),nrow=3)
-kalz(covvar)
-# example of all .1 matrix
-covvar2 <- matrix(rep(.1,9),nrow=3)
-kalz(covvar2)
-kalz <- function(covvar) {
-n <- length(covvar[1,])
-alpha <- eigen(covvar)$values[1] -10^(-2)
-kalz.ec <- solve(diag(n)-alpha*cov(covvar))%*%rep(1,n)
-return(kalz.ec)
+else stop(mess)
 }
-# example of all 1 matrix
-covvar <- matrix(rep(1,9),nrow=3)
-kalz(covvar)
-# example of all .1 matrix
-covvar2 <- matrix(rep(.1,9),nrow=3)
-kalz(covvar2)
-covvar3 <- eigen(matrix(c(1,0,0,0,1,0,0,0,1),nrow=3))
-covvar3
-covvar3 <- matrix(c(1,0,0,0,1,0,0,0,1),nrow=3)
-covvar3
-kalz(covvar3)
-covvar
-covvar2
-covvar2 <- matrix(c(1,0.1,0.1,0.1,1,0.1,0.1,0.1,1),nrow=3)
-kalz(covvar2)
-covvar2 <- matrix(c(1,0.1,0.1,0.1,1,0.1,0.1,0.1,1),nrow=3)
-kalz(covvar2)
-covvar2
-diag(covvar2) <- c(0,0,0)
-covvar2
-kalz(covvar2)
-matrix(c(1,1,-1,1,1,-1,-1,-1,1),nrow=3)
-covvar4 <- matrix(c(1,1,-1,1,1,-1,-1,-1,1),nrow=3)
-kalz(covvar4)
-covvar4
-kalz(matrix(c(1,1,-.1,1,1,-.1,-.1,-.1,1),nrow=3))
-matrix(c(1,1,-.1,1,1,-.1,-.1,-.1,1),nrow=3)
-kalz(matrix(c(1,1,.1,1,1,.1,.1,.1,1),nrow=3))
-kalz <- function(covvar) {
-n <- length(covvar[1,])
-alpha <- (eigen(covvar)$values[1])^(-1) -10^(-2)
-kalz.ec <- solve(diag(n)-alpha*cov(covvar))%*%rep(1,n)
-return(kalz.ec)
 }
-covvar <- matrix(rep(1,9),nrow=3)
-kalz(covvar)
-covvar2 <- matrix(c(1,0.1,0.1,0.1,1,0.1,0.1,0.1,1),nrow=3)
-kalz(covvar2)
-covvar3 <- matrix(c(1,0,0,0,1,0,0,0,1),nrow=3)
-kalz(covvar3)
-covvar4 <- matrix(c(1,1,-1,1,1,-1,-1,-1,1),nrow=3)
-kalz(covvar4)
-kalz(matrix(c(1,1,-.1,1,1,-.1,-.1,-.1,1),nrow=3))
-kalz(matrix(c(1,1,.1,1,1,.1,.1,.1,1),nrow=3))
-covvar4 <- matrix(c(1,1,-1,1,1,-1,-1,-1,1),nrow=3)
-kalz(covvar4)
-covvar4 <- matrix(c(1,1,-.9,1,1,-.9,-.9,-.9,1),nrow=3)
-kalz(covvar4)
-kalz(matrix(c(1,1,-.1,1,1,-.1,-.1,-.1,1),nrow=3))
-kalz(matrix(c(1,1,-1,1,1,-1,-1,-1,1),nrow=3))
-kalz(matrix(c(1,1,-1,1,1,-1,-1,-1,1),nrow=3))
-kalz(matrix(c(1,1,-.1,1,1,-.1,-.1,-.1,1),nrow=3))
-kalz(matrix(c(1,1,.1,1,1,.1,.1,.1,1),nrow=3))
-matrix(c(1,1,.1,1,1,.1,.1,.1,1),nrow=3)
-kalz(matrix(c(1,1,0,1,1,0,0,0,1),nrow=3)
-kalz(matrix(c(1,1,0,1,1,0,0,0,1),nrow=3))
-###########################################################
-kalz(matrix(c(1,1,0,1,1,0,0,0,1),nrow=3))
-matrix(c(1,1,0,1,1,0,0,0,1),nrow=3)
-###########################################################
-kalz(matrix(c(1,1,-1,1,1,-1,-1,-1,1),nrow=3))
-kalz(matrix(c(1,1,-.1,1,1,-.1,-.1,-.1,1),nrow=3))
-kalz(matrix(c(1,1,.5,1,1,.5,.5,.5,1),nrow=3))
-kalz(matrix(c(1,1,-.5,1,1,-.5,-.5,-.5,1),nrow=3))
-kalz(matrix(c(1,1,.9,1,1,.9,.9,.9,1),nrow=3))
-kalz(matrix(c(1,1,-.9,1,1,-.9,-.9,-.9,1),nrow=3))
-library(matrixcalc)
-install.packages("matrixcalc")
-is.positive.definite(matrix(c(1,1,-.1,1,1,-.1,-.1,-.1,1),nrow=3))
-is.positive.definite(matrix(c(1,1,-1,1,1,-1,-1,-1,1),nrow=3))
-library(matrixcalc)
-is.positive.definite(matrix(c(1,1,-1,1,1,-1,-1,-1,1),nrow=3))
-is.positive.definite(matrix(c(1,1,-.1,1,1,-.1,-.1,-.1,1),nrow=3))
-is.positive.definite(matrix(c(1,1,.1,1,1,.1,.1,.1,1),nrow=3))
-is.positive.definite(matrix(c(1,1,0,1,1,0,0,0,1),nrow=3))
-is.positive.definite(matrix(c(1,1,0,1,1,0,0,0,1),nrow=3))
-is.positive.definite(diag(3))
-is.positive.definite(matrix(c(1,1,-1,1,1,-1,-1,-1,1),nrow=3))
-library(mvtnorm)
-library(sna)
-library(matrixcalc)
-library(corpcor)
-chol <- cbind(c(1,0.1,0.1,0.1,0.1),c(0,0.99,0.1,0.1,0.2),c(0,0,0.98,0.4,0.5),
-c(0,0,0,0.9,0.5),c(0,0,0,0,0.7))
-covvar <- chol%*%t(chol)
-covvar
-is.positive.definite(covvar)
-eigen(covvar)
-eigen(solve(covvar))
-is.positive.definite(covvar)
-gplot(covvar,gmode="graph",edge.lwd=15,label=c(1,2,3,4,5))
-eigen(covvar)
-? gplot
-install.packages(c("JGR","Deducer","DeducerExtras"))
-library(JGR)
-JGR()
-install.packages("rJava")
-JPR()
-JGR()
-library(JGR)
-plot.lm
-library(leaps)
-library(PerformanceAnalytics)
-library(lars)
-library(robust)
-library(ellipse)
-library(MASS)
-#
-# fitMacroeconomicFactormodel
-#
-# load data from the database
-setwd("C:/Users/Yi-An Chen/Documents/R-project/factoranalytics/pkg/factorAnalytics/data")
-# data(managers.df)
-load("managers.df.rda")
-ret.assets = managers.df[,(1:6)]
-factors    = managers.df[,(7:9)]
-# fit the factor model with OLS
-setwd("C:/Users/Yi-An Chen/Documents/R-project/factoranalytics/pkg/factorAnalytics/R")
-source("fitMacroeconomicFactorModel.r")
-source("factorModelCovariance.r")
-source("factorModelSdDecomposition.r")
-source("factorModelEsDecomposition.r")
-source("factorModelVaRDecomposition.r")
-fit.macro <- fitMacroeconomicFactorModel(ret.assets,factors,fit.method="OLS", factor.set = 3,
-variable.selection="all subsets",decay.factor = 0.95)
-source("factorModelPerformanceAttribution.r")
-fm.attr <- factorModelPerformanceAttribution(fit.macro)
-fm.attr[[1]]
-fm.attr[[2]]
-fm.attr[[3]]
-fit <- fitMacroeconomicFactorModel(ret.assets,factors,fit.method="OLS", factor.set = 3,
-variable.selection="all subsets",decay.factor = 0.95)
-fm.attr <- factorModelPerformanceAttribution(fit)
-fit$ret.assets
-factors
-benchmark = managers.df[,8]
-fit$ret.assets - benchmark
-fit = fitMacroeconomicFactorModel(port.ret,fit$factors)
-port.ret =  fit$ret.assets - benchmark
-fit = fitMacroeconomicFactorModel(port.ret,fit$factors)
-fit$call
-fit <- fitMacroeconomicFactorModel(ret.assets,factors,fit.method="OLS", factor.set = 3,
-variable.selection="all subsets",decay.factor = 0.95)
-fit$call
-eval(fit$call)
-fit$call
-ret.assets =  fit$ret.assets - benchmark
-fit.1 = eval(fit$call)
-eval(fit$call)
-setwd("C:/Users/Yi-An Chen/Documents/R-project/factoranalytics/pkg/factorAnalytics/data")
-# data(managers.df)
-load("managers.df.rda")
-ret.assets = managers.df[,(1:6)]
-factors    = managers.df[,(7:9)]
-# fit the factor model with OLS
-setwd("C:/Users/Yi-An Chen/Documents/R-project/factoranalytics/pkg/factorAnalytics/R")
-source("fitMacroeconomicFactorModel.r")
-source("factorModelCovariance.r")
-source("factorModelSdDecomposition.r")
-source("factorModelEsDecomposition.r")
-source("factorModelVaRDecomposition.r")
-fit <- fitMacroeconomicFactorModel(ret.assets,factors,fit.method="OLS", factor.set = 3,
-variable.selection="all subsets",decay.factor = 0.95)
-source("factorModelPerformanceAttribution.r")
-fm.attr <- factorModelPerformanceAttribution(fit)
-fm.attr[[1]]
-source("factorModelPerformanceAttribution.r")
-fm.attr <- factorModelPerformanceAttribution(fit)
-source("factorModelPerformanceAttribution.r")
-fm.attr <- factorModelPerformanceAttribution(fit)
-source("factorModelPerformanceAttribution.r")
-fm.attr <- factorModelPerformanceAttribution(fit)
-benchmark=NULL
-if(benchmark != NULL)
+wls.classic <- function(xdf, modelterms, conlist, w) {
+assign("w", w, pos = 1)
+model <- try(lm(formula = modelterms, data = xdf, contrasts = conlist,
+weights = w, singular.ok = FALSE))
+if (is(model, "Error")) {
+mess <- geterrmessage()
+nn <- regexpr("computed fit is singular", mess)
+if (nn > 0) {
+cat("At time:", substring(mess, nn), "\n")
+model <- lm(formula = modelterms, data = xdf,
+contrasts = conlist, weights = w)
 }
-if(benchmark != NULL) {
+else stop(mess)
 }
-benchmark != NULL
-benchmark[1] != NULL
-class(benchmark)
-as.logic(benchmark)
-? logic
-? as.numeric
-as.logical(benchmark)
-(as.logical(benchmark) != NULL)
-source("factorModelPerformanceAttribution.r")
-fm.attr <- factorModelPerformanceAttribution(fit)
-source("factorModelPerformanceAttribution.r")
-fm.attr <- factorModelPerformanceAttribution(fit)
-fm.attr[[1]]
-fm.attr[[2]]
-fm.attr[[3]]
-benchmark = managers.df[,8]
-fm.attr.b <- factorModelPerformanceAttribution(fit,benchmark=benchmark)
-fm.attr.b[[1]]
-fm.attr[[1]]
-fm.attr[[2]]
-source("plot.FM.attribution.r")
-plot(fm.attr,date="2006-12-30")
-plot(fm.attr.b,date="2006-12-30")
-source("summary.FM.attribution.r")
-summary(fm.attr)
-summary(fm.attr.b)
+tstat <- rep(NA, length(model$coef))
+tstat[!is.na(model$coef)] <- summary(model, cor = FALSE)$coef[,3]
+alphaord <- order(names(model$coef))
+c(length(model$coef), model$coef[alphaord], tstat[alphaord],
+model$resid)
+}
+resids <- by(data = data, INDICES = as.numeric(data[[datevar]]),
+FUN = function(xdf, modelterms, conlist) {
+lm(formula = modelterms, data = xdf, contrasts = conlist,
+singular.ok = TRUE)$resid
+},
+modelterms = regression.formula, conlist = contrasts.list)
+resids <- apply(resids, 1, unlist)
+weights <- if (covariance == "robust")
+apply(resids, 1, scaleTau2)^2
+else apply(resids, 1, var)
+FE.hat <- by(data = data, INDICES = as.numeric(data[[datevar]]),
+FUN = wls.classic, modelterms = regression.formula,
+conlist = contrasts.list, w = weights)
+covariance = "classic"
+wls = TRUE
+regression = "classic"
+if (!wls) {
+if (regression == "robust") {
+# ols.robust
+FE.hat <- by(data = data, INDICES = as.numeric(data[[datevar]]),
+FUN = ols.robust, modelterms = regression.formula,
+conlist = contrasts.list)
+} else {
+# ols.classic
+FE.hat <- by(data = data, INDICES = as.numeric(data[[datevar]]),
+FUN = ols.classic, modelterms = regression.formula,
+conlist = contrasts.list)
+}
+} else {
+if (regression == "robust") {
+# wls.robust
+resids <- by(data = data, INDICES = as.numeric(data[[datevar]]),
+FUN = function(xdf, modelterms, conlist) {
+lmRob(modelterms, data = xdf, contrasts = conlist,
+control = lmRob.control(mxr = 200, mxf = 200,
+mxs = 200))$resid
+}, modelterms = regression.formula, conlist = contrasts.list)
+resids <- apply(resids, 1, unlist)
+weights <- if (covariance == "robust")
+apply(resids, 1, scaleTau2)^2
+else apply(resids, 1, var)
+FE.hat <- by(data = data, INDICES = as.numeric(data[[datevar]]),
+FUN = wls.robust, modelterms = regression.formula,
+conlist = contrasts.list, w = weights)
+}
+else {
+# wls.classic
+resids <- by(data = data, INDICES = as.numeric(data[[datevar]]),
+FUN = function(xdf, modelterms, conlist) {
+lm(formula = modelterms, data = xdf, contrasts = conlist,
+singular.ok = TRUE)$resid
+},
+modelterms = regression.formula, conlist = contrasts.list)
+resids <- apply(resids, 1, unlist)
+weights <- if (covariance == "robust")
+apply(resids, 1, scaleTau2)^2
+else apply(resids, 1, var)
+FE.hat <- by(data = data, INDICES = as.numeric(data[[datevar]]),
+FUN = wls.classic, modelterms = regression.formula,
+conlist = contrasts.list, w = weights)
+}
+}
+FE.hat
+FE.hat[1]
+FE.hat[[1]]
+(length(exposures.factor))
+exposures.factor
+length(levels(data[,exposures.factor]))
+(length(exposures.factor)>0)
+if (length(exposures.factor)>0) {
+numCoefs <- length(exposures.numeric) + length(levels(data[,exposures.factor]))
+ncols <- 1 + 2 * numCoefs + numAssets
+fnames <- c(exposures.numeric, paste(exposures.factor,
+levels(data[, exposures.factor]), sep = ""))
+cnames <- c("numCoefs", fnames, paste("t", fnames, sep = "."),
+assets)
+} else {
+numCoefs <- 1 + length(exposures.numeric)
+ncols <- 1 + 2 * numCoefs + numAssets
+cnames <- c("numCoefs", "(Intercept)", exposures.numeric,
+paste("t", c("(Intercept)", exposures.numeric), sep = "."),
+assets)
+}
+FE.hat.mat <- matrix(NA, ncol = ncols, nrow = numTimePoints,
+dimnames = list(as.character(timedates), cnames))
+for (i in 1:length(FE.hat)) {
+names(FE.hat[[i]])[1] <- "numCoefs"
+nc <- FE.hat[[i]][1]
+names(FE.hat[[i]])[(2 + nc):(1 + 2 * nc)] <- paste("t",
+names(FE.hat[[i]])[2:(1 + nc)], sep = ".")
+if (length(FE.hat[[i]]) != (1 + 2 * nc + numAssets))
+stop(paste("bad count in row", i, "of FE.hat"))
+names(FE.hat[[i]])[(2 + 2 * nc):(1 + 2 * nc + numAssets)] <- assets
+idx <- match(names(FE.hat[[i]]), colnames(FE.hat.mat))
+FE.hat.mat[i, idx] <- FE.hat[[i]]
+}
+coefs.names <- colnames(FE.hat.mat)[2:(1 + numCoefs)]
+# estimated factors returns ordered by time
+f.hat <- xts(x = FE.hat.mat[, 2:(1 + numCoefs)], order.by = timedates)
+# check for outlier
+gomat <- apply(coredata(f.hat), 2, function(x) abs(x - median(x,
+na.rm = TRUE)) > 4 * mad(x, na.rm = TRUE))
+if (any(gomat, na.rm = TRUE) ) {
+cat("\n\n*** Possible outliers found in the factor returns:\n\n")
+for (i in which(apply(gomat, 1, any, na.rm = TRUE))) print(f.hat[i,
+gomat[i, ], drop = FALSE])
+}
+tstats <- xts(x = FE.hat.mat[, (2 + nc):(1 + 2 * nc)], order.by = timedates)
+# residuals for every asset ordered by time
+resids <- xts(x = FE.hat.mat[, (2 + 2 * numCoefs):(1 + 2 *
+numCoefs + numAssets)], order.by = timedates)
+Cov.factors <- covClassic(coredata(f.hat), distance = FALSE,na.action = na.omit)
+resid.vars <- apply(coredata(resids), 2, var, na.rm = TRUE)
+D.hat <- if (full.resid.cov) {
+covClassic(coredata(resids), distance = FALSE, na.action = na.omit)
+}  else { diag(resid.vars) }
+B.final <- matrix(0, nrow = numAssets, ncol = numCoefs)
+colnames <-  coefs.names
+B.final
+B.final[, match("(Intercept)", colnames, 0)]
+B.final[, match("(Intercept)", colnames, 0)] <- 1
+B.final
+numeric.columns <- match(exposures.numeric, colnames, 0)
+# only take the latest beta to compute FM covariance
+B.final[, numeric.columns] <- as.matrix(data[ (data[[datevar]] == timedates[numTimePoints]), exposures.numeric])
+rownames(B.final) = assets
+colnames(B.final) = colnames(f.hat)
+B.final
+(length(exposures.factor))
+if (length(exposures.factor)>0) {
+B.final[, grep(exposures.factor, x = colnames)][cbind(seq(numAssets),
+(data[ data[[datevar]] == timedates[numTimePoints],
+exposures.factor]))] <- 1
+}
+B.final
+cov.returns <- B.final %*% Cov.factors$cov %*% t(B.final) +
+if (full.resid.cov) { D.hat$cov
+}  else { D.hat  }
+mean.cov.returns = tapply(data[[returnsvar]],data[[assetvar]], mean)
+Cov.returns <- list(cov = cov.returns, mean=mean.cov.returns, eigenvalues = eigen(cov.returns,
+only.values = TRUE, symmetric = TRUE)$values)
+# report residual covaraince if full.resid.cov is true.
+if (full.resid.cov) {
+Cov.resids <- D.hat
+}
+else {
+Cov.resids <- diag(resid.vars)
+}
+f.hat
+head(barra$factor.returns)
+barra <- fitFundamentalFactorModel(data=stacked.returns,exposure.names="SECTOR",
+datevar="DATE",returnsvar="RETURN",
+assetvar="TICKER",wls=TRUE,full.resid.cov=FALSE)
+head(barra$factor.returns)
+head(f.hat)
+(!(length(exposures.factor)>0))
+setwd("C:/Users/Yi-An Chen/Documents/R-project/returnanalytics/pkg/FactorAnalytics/R")
+source(fitFundamentalFactorModel)
+source("fitFundamentalFactorModel.r")
+barra <- fitFundamentalFactorModel(data=stacked.returns,exposure.names="SECTOR",
+datevar="DATE",returnsvar="RETURN",
+assetvar="TICKER",wls=TRUE,full.resid.cov=FALSE)
+names(barra)
+head(barra$factor.returns)
+barra$beta
+data(berndtInvest)
+ts.berndt<-xts(berndtInvest[,-1],as.Date(berndtInvest[,1]))
+berndt<-ts.berndt['1978/1987']
+returns<-berndt[,c(-10,-17)]
+tickers <- names(returns)
+num.tickers <- length(tickers)
+dates <- index(returns)
+num.dates <- length(dates)
+sector<-c("OTHER","OTHER","OIL","TECH","TECH","OIL","OTHER","OTHER",
+"TECH","OIL","OIL","OTHER","TECH","OIL","OTHER")
+stacked.returns <- data.frame(
+DATE=rep(dates,num.tickers),
+TICKER=rep(tickers,each=num.dates),
+RETURN=c(coredata(returns)),
+SECTOR=rep(sector,each=num.dates),
+stringsAsFactors=FALSE)
+head(stacked.returns)
+setwd("C:/Users/Yi-An Chen/Documents/R-project/returnanalytics/pkg/FactorAnalytics/R")
+barra <- fitFundamentalFactorModel(data=stacked.returns,exposure.names="SECTOR",
+datevar="DATE",returnsvar="RETURN",
+assetvar="TICKER",wls=TRUE,full.resid.cov=FALSE)
+names(barra)
+head(barra$factor.returns)
+setwd("C:/Users/Yi-An Chen/Documents/R-project/returnanalytics/pkg/FactorAnalytics/R")
+source("fitFundamentalFactorModel.r")
+barra <- fitFundamentalFactorModel(data=stacked.returns,exposure.names="SECTOR",
+datevar="DATE",returnsvar="RETURN",
+assetvar="TICKER",wls=TRUE,full.resid.cov=FALSE)
+names(barra)
+head(barra$factor.returns)
+barra$beta
+barra.cov <- factorModelCovariance(barra$beta, barra$factor.cov$cov, barra$resid.variance)
+barra.cor <- cov2cor(barra.cov)
+round(barra.cor,2)
+returns = berndtInvest[,-c(1,11,18)]
+n.stocks = ncol(returns)
+tech.dum = oil.dum = other.dum = matrix(0,n.stocks,1)
+tech.dum[c(4,5,9,13),] = 1
+oil.dum[c(3,6,10,11,14),] = 1
+other.dum = 1 - tech.dum - oil.dum
+B = cbind(tech.dum,oil.dum,other.dum)
+dimnames(B) = list(colnames(returns),c("TECH","OIL","OTHER"))
+B
+returns = t(returns)
+barra_ols = lm(returns ~ -1 + B)
+sbarra_ols = summary(barra_ols)
+e.hat = resid(barra_ols)
+e.var = apply(e.hat,1,var)
+e.var.inv = e.var^-1
+barra_fgls = lm(returns ~ -1 + B, weights = e.var.inv)
+sbarra_fgls = summary(barra_fgls)
+f.hat = coef(barra_fgls)
+f.hat
+e.hat_fgls = resid(barra_fgls)
+e.var_fgls = apply(e.hat_fgls,1,var)
+cov.ind.lm = B %*% var(t(f.hat)) %*% t(B) + diag(e.var_fgls)
+all.equal(cov.ind,cov.ind.lm)
+cor.ind.lm <- cov2cor(cov.ind.lm)
+round(cor.ind.lm,2)
+cov.ind.lm = B %*% var(t(f.hat)) %*% t(B) + diag(e.var_fgls)
+all.equal(cov.ind,cov.ind.lm)
+F.hat = solve(crossprod(B))%*%t(B)%*%returns
+# compute N x T matrix of industry factor model residuals
+E.hat = returns - B%*%F.hat
+# compute residual variances from time series of errors
+diagD.hat = apply(E.hat,1,var)
+Dinv.hat = diag(diagD.hat^(-1))
+# multivariate FGLS regression to estimate K x T matrix of factor returns
+H = solve(t(B)%*%Dinv.hat%*%B)%*%t(B)%*%Dinv.hat
+round(H[,1:8],5)
+# note: rows of H sum to one
+apply(H,1,sum)
+# create factor mimicking portfolios
+F.hat = H%*%returns
+E.hat = returns - B%*%F.hat
+diagD.hat = apply(E.hat,1,var)
+F.hat = t(F.hat)
+all.equal(cov.ind,cov.ind.lm)
+cov.ind = B%*%var(F.hat)%*%t(B) + diag(diagD.hat)
+sd = sqrt(diag(cov.ind))
+cor.ind = cov.ind/outer(sd,sd)
+cor.samp <- cor(t(returns))
+all.equal(cov.ind,cov.ind.lm)
+cor.ind.lm <- cov2cor(cov.ind.lm)
+round(cor.ind.lm,2)
+barra_ols = lm(returns ~ -1 + B)
+sbarra_ols = summary(barra_ols)
+e.hat = resid(barra_ols)
+e.var = apply(e.hat,1,var)
+# e.var.inv = e.var^-1
+e.var.inv = e.var
+barra_fgls = lm(returns ~ -1 + B, weights = e.var.inv)
+sbarra_fgls = summary(barra_fgls)
+f.hat = coef(barra_fgls)
+e.hat_fgls = resid(barra_fgls)
+e.var_fgls = apply(e.hat_fgls,1,var)
+cov.ind.lm = B %*% var(t(f.hat)) %*% t(B) + diag(e.var_fgls)
+all.equal(cov.ind,cov.ind.lm)
+barra.cov <- factorModelCovariance(barra$beta, barra$factor.cov$cov, barra$resid.variance)
+all.equal(cov.ind.lm,barra.cov)
+? lm
+setwd("C:/Users/Yi-An Chen/Documents/R-project/returnanalytics/pkg/FactorAnalytics/R")
+source("fitFundamentalFactorModel.r")
+barra <- fitFundamentalFactorModel(data=stacked.returns,exposure.names="SECTOR",
+datevar="DATE",returnsvar="RETURN",
+assetvar="TICKER",wls=TRUE,full.resid.cov=FALSE)
+head(barra$factor.returns)
+barra.cov <- factorModelCovariance(barra$beta, barra$factor.cov$cov, barra$resid.variance)
+barra.cor <- cov2cor(barra.cov)
+round(barra.cor,2)
+barra_ols = lm(returns ~ -1 + B)
+sbarra_ols = summary(barra_ols)
+e.hat = resid(barra_ols)
+e.var = apply(e.hat,1,var)
+e.var.inv = e.var^-1
+barra_fgls = lm(returns ~ -1 + B, weights = e.var.inv)
+sbarra_fgls = summary(barra_fgls)
[TRUNCATED]

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