[Returnanalytics-commits] r2200 - pkg/FactorAnalytics/R

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Mon Jul 23 18:37:26 CEST 2012


Author: chenyian
Date: 2012-07-23 18:37:26 +0200 (Mon, 23 Jul 2012)
New Revision: 2200

Added:
   pkg/FactorAnalytics/R/plot.FundamentalFactorModel.r
   pkg/FactorAnalytics/R/plot.StatFactorModel.r
   pkg/FactorAnalytics/R/print.StatFactorModel.r
Modified:
   pkg/FactorAnalytics/R/fitFundamentalFactorModel.R
   pkg/FactorAnalytics/R/fitMacroeconomicFactorModel.R
   pkg/FactorAnalytics/R/fitStatisticalFactorModel.R
   pkg/FactorAnalytics/R/plot.MacroFactorModel.r
Log:


Modified: pkg/FactorAnalytics/R/fitFundamentalFactorModel.R
===================================================================
--- pkg/FactorAnalytics/R/fitFundamentalFactorModel.R	2012-07-23 06:18:37 UTC (rev 2199)
+++ pkg/FactorAnalytics/R/fitFundamentalFactorModel.R	2012-07-23 16:37:26 UTC (rev 2200)
@@ -4,7 +4,8 @@
           datevar = "DATE", assetvar = "PERMNO", returnsvar = "RETURN", 
           tickersvar = "TICKER.x") {
   
-
+require(zoo)
+require(robust)
 ## input
 ##  
 ## fulldata               : data.frame. data stacked by dates
@@ -32,15 +33,15 @@
 ## 
 ## cov.returns            : covariance information for asset returns, includes
 ##                          covariance, mean, eigenvalus
-## cov.factor.rets            : covariance information for factor returns, includes
+## cov.factor.rets        : covariance information for factor returns, includes
 ##                          covariance and mean
 ## cov.resids             : covariance information for model residuals, includes
 ##                          covariance and mean
 ## resid.vars             : list of information regarding model residuals variance
-## factor.rets                : zoo time series object of factor returns
+## factor.rets            : zoo time series object of factor returns
 ## resids                 : zoo time series object of model residuals
 ## tstats                 : zoo time series object of model t-statistics
-##returns.data            : data.frame of return data including RETURN, DATE,PERMNO
+## returns.data            : data.frame of return data including RETURN, DATE,PERMNO
 ## exposure.data          : data.frame of exposure data including DATE, PERMNO
 ## assets                 : character vector of PERMNOs used in the model
 ## tickers                : character vector of tickers used in the model
@@ -338,6 +339,7 @@
                    assets = assets, 
                    tickers = tickers, 
                    call = this.call)
+    class(output) <- "FundamentalFactorModel"
     return(output)
 }
 

Modified: pkg/FactorAnalytics/R/fitMacroeconomicFactorModel.R
===================================================================
--- pkg/FactorAnalytics/R/fitMacroeconomicFactorModel.R	2012-07-23 06:18:37 UTC (rev 2199)
+++ pkg/FactorAnalytics/R/fitMacroeconomicFactorModel.R	2012-07-23 16:37:26 UTC (rev 2200)
@@ -38,7 +38,11 @@
 ## r2.vec              vector of R-square values
 ## residVars.vec       vector of residual variances
 ## call               function call.  
-
+  require(leaps)
+  require(lars)
+  require(robust)
+  require(ellipse)
+  require(MASS)
 this.call <- match.call()
   
 if (is.data.frame(ret.assets) & is.data.frame(factors) ) {

Modified: pkg/FactorAnalytics/R/fitStatisticalFactorModel.R
===================================================================
--- pkg/FactorAnalytics/R/fitStatisticalFactorModel.R	2012-07-23 06:18:37 UTC (rev 2199)
+++ pkg/FactorAnalytics/R/fitStatisticalFactorModel.R	2012-07-23 16:37:26 UTC (rev 2200)
@@ -55,6 +55,7 @@
 		stop("Invalid choice for optional argument method.")
 	}
  return(ans)
+    
 }
 
  
@@ -167,9 +168,21 @@
 	dimnames(f) <- list(dimnames(x)[[1]], paste("F", 1:k, sep = "."))
 	dimnames(Omega) <- list(x.names, x.names)
 	names(alpha) <- x.names
+  # create lm list for plot
+  reg.list = list()
+  for (i in x.names) {
+    reg.df = as.data.frame(cbind(x[,i],f))
+    colnames(reg.df)[1] <- i
+    fm.formula = as.formula(paste(i,"~", ".", sep=" "))
+    fm.fit = lm(fm.formula, data=reg.df)
+    reg.list[[i]] = fm.fit
+    }
 	ans <-  list(factors = f, loadings = B, k = k, alpha = alpha, Omega = Omega,
-	            	r2 = r2, eigen = eigen.tmp$values, residuals=tmp)
- return(ans)
+	            	r2 = r2, eigen = eigen.tmp$values, residuals=tmp, asset.ret = x,
+               asset.fit=reg.list)
+ 
+  return(ans)
+ 
 }
 
   # funciont of apca
@@ -217,7 +230,7 @@
 	res <- t(t(x) - alpha) - f %*% B
 	r2 <- (1 - colSums(res^2)/colSums(xc^2))
   ans <- 	list(factors = f, loadings = B, k = k, alpha = alpha, Omega = Omega,
-		           r2 = r2, eigen = eig.tmp$values, residuals=res)
+		           r2 = r2, eigen = eig.tmp$values, residuals=res,asset.ret = x)
  return(ans)
 }
 
@@ -265,12 +278,14 @@
 	}	else {
 		mimic <- qr.solve(x, f)
 	}
-	mimic <- t(t(mimic)/colSums(mimic))
+	
+  mimic <- t(t(mimic)/colSums(mimic))
 	dimnames(mimic)[[1]] <- dimnames(x)[[2]]
+  
   ans$mimic <- mimic
   ans$residVars.vec <- apply(ans$residuals,2,var)
   ans$call <- call
-  oldClass(ans) <- "sfm"
+class(ans) <- "StatFactorModel"
   return(ans)
 }
 

Added: pkg/FactorAnalytics/R/plot.FundamentalFactorModel.r
===================================================================
--- pkg/FactorAnalytics/R/plot.FundamentalFactorModel.r	                        (rev 0)
+++ pkg/FactorAnalytics/R/plot.FundamentalFactorModel.r	2012-07-23 16:37:26 UTC (rev 2200)
@@ -0,0 +1,58 @@
+# plot.FundamentalFactorModel.r
+# Yi-An Chen
+# 7/16/2012
+
+plot.FundamentalFactorModel <- 
+function(fund.fit,which.plot=c("none","1L","2L","3L","4L","5L","6L","7L"),max.show=12,
+         plot.single=FALSE, fundId, fundName="TBA",
+         which.plot.single=c("none","1L","2L","3L","4L","5L","6L","7L","8L",
+                             "9L","10L","11L","12L","13L") ) {
+require(ellipse)
+  
+  if (plot.single==TRUE) {
+  
+  } else {    
+    which.plot<-which.plot[1]
+    
+    if(which.plot=='none') 
+      which.plot<-menu(c("Factor returns",
+                         "Residual plots",
+                         "Variance of Residuals",
+                         "Factor Model Correlation"),
+                       title="Factor Analytics Plot \nMake a plot selection (or 0 to exit):\n") 
+    
+    
+    n <- length(fund.fit$asset)
+    if (n >= max.show) {
+      cat(paste("numbers of assets are greater than",max.show,", show only first",
+                max.show,"assets",sep=" "))
+      n <- max.show 
+    }
+    switch(which.plot,
+           
+           "1L" = {
+            plot(fund.fit$factor.rets,main="Factor Returns") 
+             
+           }, 
+          "2L" ={
+            plot(fund.fit$resids[,c(1:n)],main="Residuals")
+           },
+           "3L" = {
+             barplot(fund.fit$resid.vars[c(1:n)])  
+           },    
+           
+           "4L" = {
+             cor.fm = cov2cor(fund.fit$cov.returns$cov)
+             rownames(cor.fm) = colnames(cor.fm)
+             ord <- order(cor.fm[1,])
+             ordered.cor.fm <- cor.fm[ord, ord]
+             plotcorr(ordered.cor.fm[c(1:n),c(1:n)], col=cm.colors(11)[5*ordered.cor.fm + 6])
+           },
+           
+           invisible()       
+    )         
+  }
+  
+  
+} 
+  
\ No newline at end of file

Modified: pkg/FactorAnalytics/R/plot.MacroFactorModel.r
===================================================================
--- pkg/FactorAnalytics/R/plot.MacroFactorModel.r	2012-07-23 06:18:37 UTC (rev 2199)
+++ pkg/FactorAnalytics/R/plot.MacroFactorModel.r	2012-07-23 16:37:26 UTC (rev 2200)
@@ -1,25 +1,200 @@
-plot.MacroFactorModel <- 
-  function(fit.macro,colorset=c(1:12),legend.loc=NULL,which=c("none","1L","2L","3L",
-                                                              "4L","5L")) {
-    which<-which[1]
+ plot.MacroFactorModel <- 
+  function(fit.macro,colorset=c(1:12),legend.loc=NULL,
+           which.plot=c("none","1L","2L","3L","4L","5L","6L","7L"),max.show=6,
+           plot.single=FALSE, fundId, fundName="TBA",which.plot.single=c("none","1L","2L","3L","4L","5L","6L",
+                                                                  "7L","8L","9L","10L","11L","12L","13L")) {
+      require(zoo)
+      require(PerformanceAnalytics)
+      require(strucchange)
     
-    if(which=='none') 
-      which<-menu(c("Fitted factor returns","FM Correlation",
-                    "Factor Contributions to SD",
-                    "Factor Contributions to ES",
-                    "Factor Contributions to VaR"),
-                  title="Factor Analytics Plot") 
+    if (plot.single==TRUE) {
+      ## inputs:
+      ## fit.macro        lm object summarizing factor model fit. It is assumed that
+      ##                  time series date information is included in the names component
+      ##                  of the residuals, fitted and model components of the object.
+      ## fundId           charater. The name of the single asset to be ploted.            
+      ## fundName         TBA
+      ## which.plot.single       integer indicating which plot to create:
+      ##                  1     time series plot of actual and fitted values
+      ##                  2     time series plot of residuals with standard error bands
+      ##                  3     time series plot of squared residuals
+      ##                  4     time series plot of absolute residuals
+      ##                  5     SACF and PACF of residuals
+      ##                  6     SACF and PACF of squared residuals
+      ##                  7     SACF and PACF of absolute residuals
+      ##                  8     histogram of residuals with normal curve overlayed
+      ##                  9     normal qq-plot of residuals
+      ##                  10    CUSUM plot of recursive residuals
+      ##                  11    CUSUM plot of OLS residuals
+      ##                  12    CUSUM plot of recursive estimates relative to full sample estimates
+      ##                  13    rolling estimates over 24 month window
+      which.plot.single<-which.plot.single[1]
+      fit.lm = fit.macro$asset.fit[[fundId]]
+      
+      if (!(class(fit.lm) == "lm"))
+        stop("Must pass a valid lm object")
+      
+      ## extract information from lm object
+        
+      factorNames = colnames(fit.lm$model)[-1]
+      fit.formula = as.formula(paste(fundId,"~", paste(factorNames, collapse="+"), sep=" "))
+      residuals.z = zoo(residuals(fit.lm), as.Date(names(residuals(fit.lm))))
+      fitted.z = zoo(fitted(fit.lm), as.Date(names(fitted(fit.lm))))
+      actual.z = zoo(fit.lm$model[,1], as.Date(rownames(fit.lm$model)))
+      tmp.summary = summary(fit.lm)
+       
+        
+      if (which.plot.single=="none")
+      which.plot.single<-menu(c("time series plot of actual and fitted values",
+                         "time series plot of residuals with standard error bands",
+                         "time series plot of squared residuals",
+                         "time series plot of absolute residuals",
+                         "SACF and PACF of residuals",
+                         "SACF and PACF of squared residuals",
+                         "SACF and PACF of absolute residuals",
+                         "histogram of residuals with normal curve overlayed",
+                         "normal qq-plot of residuals",
+                         "CUSUM plot of recursive residuals",
+                         "CUSUM plot of OLS residuals",
+                         "CUSUM plot of recursive estimates relative to full sample estimates",
+                         "rolling estimates over 24 month window"),
+                       title="\nMake a plot selection (or 0 to exit):\n")
+      switch(which.plot.single,
+             "1L" =  {
+        ##  time series plot of actual and fitted values
+        plot(actual.z, main=fundName, ylab="Monthly performance", lwd=2, col="black")
+        lines(fitted.z, lwd=2, col="blue")
+        abline(h=0)
+        legend(x="bottomleft", legend=c("Actual", "Fitted"), lwd=2, col=c("black","blue"))
+      }, 
+      
+             "2L" = {
+        ## time series plot of residuals with standard error bands
+        plot(residuals.z, main=fundName, ylab="Monthly performance", lwd=2, col="black")
+        abline(h=0)
+        abline(h=2*tmp.summary$sigma, lwd=2, lty="dotted", col="red")
+        abline(h=-2*tmp.summary$sigma, lwd=2, lty="dotted", col="red")
+        legend(x="bottomleft", legend=c("Residual", "+/ 2*SE"), lwd=2,
+               lty=c("solid","dotted"), col=c("black","red"))
+      },
+             "3L" = {
+        ## time series plot of squared residuals
+        plot(residuals.z^2, main=fundName, ylab="Squared residual", lwd=2, col="black")
+        abline(h=0)
+        legend(x="topleft", legend="Squared Residuals", lwd=2, col="black")
+      },
+             "4L" = {
+        ## time series plot of absolute residuals
+        plot(abs(residuals.z), main=fundName, ylab="Absolute residual", lwd=2, col="black")
+        abline(h=0)
+        legend(x="topleft", legend="Absolute Residuals", lwd=2, col="black")
+      },
+             "5L" = {
+        ## SACF and PACF of residuals
+        chart.ACFplus(residuals.z, main=paste("Residuals: ", fundName, sep=""))
+      },
+             "6L" = {
+        ## SACF and PACF of squared residuals
+        chart.ACFplus(residuals.z^2, main=paste("Residuals^2: ", fundName, sep=""))
+      },
+             "7L" = {
+        ## SACF and PACF of absolute residuals
+        chart.ACFplus(abs(residuals.z), main=paste("|Residuals|: ", fundName, sep=""))
+      },
+             "8L" = {
+        ## histogram of residuals with normal curve overlayed
+        chart.Histogram(residuals.z, methods="add.normal", main=paste("Residuals: ", fundName, sep=""))
+      },
+             "9L" = {
+        ##  normal qq-plot of residuals
+        chart.QQPlot(residuals.z, envelope=0.95, main=paste("Residuals: ", fundName, sep=""))
+      },
+             "10L"= {
+        ##  CUSUM plot of recursive residuals
+   if (as.character(fit.macro$call["fit.method"]) == "OLS") {
+        cusum.rec = efp(fit.formula, type="Rec-CUSUM", data=fit.lm$model)
+        plot(cusum.rec, sub=fundName)
+   } else 
+     stop("CUMSUM applies only on OLS method")
+      },
+             "11L"= {
+        ##  CUSUM plot of OLS residuals
+               if (as.character(fit.macro$call["fit.method"]) == "OLS") {        
+        cusum.ols = efp(fit.formula, type="OLS-CUSUM", data=fit.lm$model)
+        plot(cusum.ols, sub=fundName)
+               } else 
+                 stop("CUMSUM applies only on OLS method")   
+      },
+             "12L"= {
+        ##  CUSUM plot of recursive estimates relative to full sample estimates
+               if (as.character(fit.macro$call["fit.method"]) == "OLS") {        
+        cusum.est = efp(fit.formula, type="fluctuation", data=fit.lm$model)
+        plot(cusum.est, functional=NULL, sub=fundName)
+               } else 
+                 stop("CUMSUM applies only on OLS method")
+      },
+             "13L"= {
+        ##  rolling regression over 24 month window
+    if (as.character(fit.macro$call["fit.method"]) == "OLS") {   
+          rollReg <- function(data.z, formula) {
+          coef(lm(formula, data = as.data.frame(data.z)))  
+        }
+        reg.z = zoo(fit.lm$model, as.Date(rownames(fit.lm$model)))
+        rollReg.z = rollapply(reg.z, FUN=rollReg, fit.formula, width=24, by.column = FALSE, 
+                              align="right")
+        plot(rollReg.z, main=paste("24-month rolling regression estimates:", fundName, sep=" "))
+    } else if (as.character(fit.macro$call["fit.method"]) == "DLS") {
+      decay.factor <- as.numeric(as.character(fit.macro$call["decay.factor"]))
+      t.length <- 24
+      w <- rep(decay.factor^(t.length-1),t.length)
+      for (k in 2:t.length) {
+        w[k] = w[k-1]/decay.factor 
+      }   
+      w <- w/sum(w)
+      rollReg <- function(data.z, formula,w) {
+        coef(lm(formula,weight=w, data = as.data.frame(data.z)))  
+      }
+      reg.z = zoo(fit.lm$model[-length(fit.lm$model)], as.Date(rownames(fit.lm$model)))
+      factorNames = colnames(fit.lm$model)[c(-1,-length(fit.lm$model))]
+      fit.formula = as.formula(paste(fundId,"~", paste(factorNames, collapse="+"), sep=" "))
+      rollReg.z = rollapply(reg.z, FUN=rollReg, fit.formula,w, width=24, by.column = FALSE, 
+                            align="right")
+      plot(rollReg.z, main=paste("24-month rolling regression estimates:", fundName, sep=" ")) 
+    } 
+        },
+             invisible()
+             )
+      
+      
+      
+    } else {    
+    which.plot<-which.plot[1]
     
+    if(which.plot=='none') 
+      which.plot<-menu(c("Fitted factor returns",
+                         "R square",
+                         "Variance of Residuals",
+                         "FM Correlation",
+                         "Factor Contributions to SD",
+                        "Factor Contributions to ES",
+                         "Factor Contributions to VaR"),
+                  title="Factor Analytics Plot \nMake a plot selection (or 0 to exit):\n") 
+    
     variable.selection = fit.macro$variable.selection
     manager.names = colnames(fit.macro$ret.assets)
     factor.names  = colnames(fit.macro$factors)
-    managers.df   = cbind(ret.assets,factors)
+    managers.df   = cbind(fit.macro$ret.assets,fit.macro$factors)
     cov.factors = var(fit.macro$factors)
     n <- length(manager.names)
     
-    
-    switch(which,"1L" = {
-          
+    switch(which.plot,
+           
+           "1L" = {
+     if (n >= max.show) {
+      cat(paste("numbers of assets are greater than",max.show,", show only first",
+                max.show,"assets",sep=" "))
+    n <- max.show 
+     }      
     par(mfrow=c(n/2,2))
     if (variable.selection == "lar" || variable.selection == "lasso") {
      for (i in 1:n) {
@@ -41,8 +216,14 @@
     }
     par(mfrow=c(1,1))
     },
+      "2L" ={
+      barplot(fit.macro$r2.vec)
+     },
+      "3L" = {
+      barplot(fit.macro$residVars.vec)  
+      },    
            
-     "2L" = {
+     "4L" = {
       cov.fm<- factorModelCovariance(fit.macro$beta.mat,var(fit.macro$factors),fit.macro$residVars.vec)    
       cor.fm = cov2cor(cov.fm)
       rownames(cor.fm) = colnames(cor.fm)
@@ -50,7 +231,7 @@
       ordered.cor.fm <- cor.fm[ord, ord]
       plotcorr(ordered.cor.fm, col=cm.colors(11)[5*ordered.cor.fm + 6])
            },
-    "3L" = {
+    "5L" = {
        factor.sd.decomp.list = list()
        for (i in manager.names) {
          factor.sd.decomp.list[[i]] =
@@ -70,9 +251,27 @@
                  col=c(1:50) )
       
     },
-     "4L"={
-       
-       factor.es.decomp.list = list()
+     "6L"={
+        factor.es.decomp.list = list()
+       if (variable.selection == "lar" || variable.selection == "lasso") {
+        
+         for (i in manager.names) {
+           idx = which(!is.na(managers.df[,i]))
+           alpha = fit.macro$alpha.vec[i]
+           beta = as.matrix(fit.macro$beta.mat[i,])        
+           fitted = alpha+as.matrix(fit.macro$factors)%*%beta
+           residual = fit.macro$ret.assets[,i]-fitted
+           tmpData = cbind(managers.df[idx,i], managers.df[idx,factor.names],
+                           (residual[idx,]/sqrt(fit.macro$residVars.vec[i])) )
+           colnames(tmpData)[c(1,length(tmpData))] = c(i, "residual")
+           factor.es.decomp.list[[i]] = 
+          factorModelEsDecomposition(tmpData, 
+                                        fit.macro$beta.mat[i,],
+                                        fit.macro$residVars.vec[i], tail.prob=0.05)
+           
+         }
+       } else {
+            
        for (i in manager.names) {
          # check for missing values in fund data
          idx = which(!is.na(managers.df[,i]))
@@ -84,7 +283,7 @@
                                       fit.macro$beta.mat[i,],
                                       fit.macro$residVars.vec[i], tail.prob=0.05)
        }
-             
+       }     
        
        # stacked bar charts of percent contributions to SD
        getCETL = function(x) {
@@ -97,9 +296,28 @@
                legend.text=T, args.legend=list(x="topleft"),
                col=c(1:50) ) 
      },
-    "5L" ={
+    "7L" ={
       
       factor.VaR.decomp.list = list()
+      
+      if (variable.selection == "lar" || variable.selection == "lasso") {
+        
+        for (i in manager.names) {
+          idx = which(!is.na(managers.df[,i]))
+          alpha = fit.macro$alpha.vec[i]
+          beta = as.matrix(fit.macro$beta.mat[i,])        
+          fitted = alpha+as.matrix(fit.macro$factors)%*%beta
+          residual = fit.macro$ret.assets[,i]-fitted
+          tmpData = cbind(managers.df[idx,i], managers.df[idx,factor.names],
+                          (residual[idx,]/sqrt(fit.macro$residVars.vec[i])) )
+          colnames(tmpData)[c(1,length(tmpData))] = c(i, "residual")
+          factor.VaR.decomp.list[[i]] = 
+            factorModelVaRDecomposition(tmpData, 
+                                       fit.macro$beta.mat[i,],
+                                       fit.macro$residVars.vec[i], tail.prob=0.05)
+          
+        }
+      } else {
       for (i in manager.names) {
         # check for missing values in fund data
         idx = which(!is.na(managers.df[,i]))
@@ -112,8 +330,8 @@
                                      fit.macro$residVars.vec[i], tail.prob=0.05,
                                       VaR.method="HS")
       }
+      }
       
-      
       # stacked bar charts of percent contributions to SD
       getCVaR = function(x) {
         x$cVaR.fm
@@ -124,8 +342,9 @@
       barplot(cr.VaR, main="Factor Contributions to VaR",
               legend.text=T, args.legend=list(x="topleft"),
               col=c(1:50) )
+    },
+    invisible()       
+    )         
     }       
-    )         
-           
  
 }
\ No newline at end of file

Added: pkg/FactorAnalytics/R/plot.StatFactorModel.r
===================================================================
--- pkg/FactorAnalytics/R/plot.StatFactorModel.r	                        (rev 0)
+++ pkg/FactorAnalytics/R/plot.StatFactorModel.r	2012-07-23 16:37:26 UTC (rev 2200)
@@ -0,0 +1,323 @@
+plot.StatFactorModel <-
+function(x, variables, cumulative = TRUE, style = "bar",
+         which.plot = c("none","1L","2L","3L","4L","5L","6L","7L","8L"),
+         hgrid = FALSE, vgrid = FALSE,plot.single=FALSE, fundId, fundName="TBA",
+         which.plot.single=c("none","1L","2L","3L","4L","5L","6L",
+                             "7L","8L","9L","10L","11L","12L","13L"), ...)
+{
+  require(strucchange)
+  #
+  # beginning of funciton screenplot
+  #
+  screeplot<-
+  function(mf, variables, cumulative = TRUE, style = "bar", main = "", ...)
+  {
+    vars <- mf$eigen
+    n90 <- which(cumsum(vars)/sum(vars) > 0.9)[1]
+    if(missing(variables)) {
+      variables <- 1:max(mf$k, min(10, n90))
+    }
+    istyle <- charmatch(style, c("bar", "lines"), nomatch = NA)
+    if(is.na(istyle) || istyle <= 1)
+      style <- "bar"
+    else {
+      style <- "lines"
+    }
+    if(style == "bar") {
+      loc <- barplot(vars[variables], names = paste("F", variables,
+                                                    sep = "."), main = main, ylab = "Variances", ...)
+    }
+    else {
+      loc <- 1:length(variables)
+      plot(loc, vars[variables], type = "b", axes = F, main = main,
+           ylab = "Variances", xlab = "")
+      axis(2)
+      axis(1, at = loc, labels = paste("F", variables, sep = "."))
+    }
+    if(cumulative) {
+      cumv <- (cumsum(vars)/sum(vars))[variables]
+      text(loc, vars[variables] + par("cxy")[2], as.character(signif(
+        cumv, 3)))
+    }
+    invisible(loc)
+  }
+  #
+  # end of screenplot
+  #
+  
+  if (plot.single==TRUE) {
+    ## inputs:
+    ## x               lm object summarizing factor model fit. It is assumed that
+    ##                  time series date information is included in the names component
+    ##                  of the residuals, fitted and model components of the object.
+    ## fundId           charater. The name of the single asset to be ploted.            
+    ## fundName         TBA
+    ## which.plot.single       integer indicating which plot to create:
+    ##                  1     time series plot of actual and fitted values
+    ##                  2     time series plot of residuals with standard error bands
+    ##                  3     time series plot of squared residuals
+    ##                  4     time series plot of absolute residuals
+    ##                  5     SACF and PACF of residuals
+    ##                  6     SACF and PACF of squared residuals
+    ##                  7     SACF and PACF of absolute residuals
+    ##                  8     histogram of residuals with normal curve overlayed
+    ##                  9     normal qq-plot of residuals
+    ##                  10    CUSUM plot of recursive residuals
+    ##                  11    CUSUM plot of OLS residuals
+    ##                  12    CUSUM plot of recursive estimates relative to full sample estimates
+    ##                  13    rolling estimates over 24 month window
+    which.plot.single<-which.plot.single[1]
+    fit.lm = x$asset.fit[[fundId]]
+    
+    if (!(class(fit.lm) == "lm"))
+      stop("Must pass a valid lm object")
+    
+    ## extract information from lm object
+    
+    factorNames = colnames(fit.lm$model)[-1]
+    fit.formula = as.formula(paste(fundId,"~", paste(factorNames, collapse="+"), sep=" "))
+    residuals.z = zoo(residuals(fit.lm), as.Date(names(residuals(fit.lm))))
+    fitted.z = zoo(fitted(fit.lm), as.Date(names(fitted(fit.lm))))
+    actual.z = zoo(fit.lm$model[,1], as.Date(rownames(fit.lm$model)))
+    tmp.summary = summary(fit.lm)
+    
+    
+    if (which.plot.single=="none")
+      which.plot.single<-menu(c("time series plot of actual and fitted values",
+                                "time series plot of residuals with standard error bands",
+                                "time series plot of squared residuals",
+                                "time series plot of absolute residuals",
+                                "SACF and PACF of residuals",
+                                "SACF and PACF of squared residuals",
+                                "SACF and PACF of absolute residuals",
+                                "histogram of residuals with normal curve overlayed",
+                                "normal qq-plot of residuals",
+                                "CUSUM plot of recursive residuals",
+                                "CUSUM plot of OLS residuals",
+                                "CUSUM plot of recursive estimates relative to full sample estimates",
+                                "rolling estimates over 24 month window"),
+                              title="\nMake a plot selection (or 0 to exit):\n")
+    switch(which.plot.single,
+           "1L" =  {
+             ##  time series plot of actual and fitted values
+             plot(actual.z, main=fundName, ylab="Monthly performance", lwd=2, col="black")
+             lines(fitted.z, lwd=2, col="blue")
+             abline(h=0)
+             legend(x="bottomleft", legend=c("Actual", "Fitted"), lwd=2, col=c("black","blue"))
+           }, 
+           
+           "2L" = {
+             ## time series plot of residuals with standard error bands
+             plot(residuals.z, main=fundName, ylab="Monthly performance", lwd=2, col="black")
+             abline(h=0)
+             abline(h=2*tmp.summary$sigma, lwd=2, lty="dotted", col="red")
+             abline(h=-2*tmp.summary$sigma, lwd=2, lty="dotted", col="red")
+             legend(x="bottomleft", legend=c("Residual", "+/ 2*SE"), lwd=2,
+                    lty=c("solid","dotted"), col=c("black","red"))
+           },
+           "3L" = {
+             ## time series plot of squared residuals
+             plot(residuals.z^2, main=fundName, ylab="Squared residual", lwd=2, col="black")
+             abline(h=0)
+             legend(x="topleft", legend="Squared Residuals", lwd=2, col="black")
+           },
+           "4L" = {
+             ## time series plot of absolute residuals
+             plot(abs(residuals.z), main=fundName, ylab="Absolute residual", lwd=2, col="black")
+             abline(h=0)
+             legend(x="topleft", legend="Absolute Residuals", lwd=2, col="black")
+           },
+           "5L" = {
+             ## SACF and PACF of residuals
+             chart.ACFplus(residuals.z, main=paste("Residuals: ", fundName, sep=""))
+           },
+           "6L" = {
+             ## SACF and PACF of squared residuals
+             chart.ACFplus(residuals.z^2, main=paste("Residuals^2: ", fundName, sep=""))
+           },
+           "7L" = {
+             ## SACF and PACF of absolute residuals
+             chart.ACFplus(abs(residuals.z), main=paste("|Residuals|: ", fundName, sep=""))
+           },
+           "8L" = {
+             ## histogram of residuals with normal curve overlayed
+             chart.Histogram(residuals.z, methods="add.normal", main=paste("Residuals: ", fundName, sep=""))
+           },
+           "9L" = {
+             ##  normal qq-plot of residuals
+             chart.QQPlot(residuals.z, envelope=0.95, main=paste("Residuals: ", fundName, sep=""))
+           },
+           "10L"= {
+             ##  CUSUM plot of recursive residuals
+            
+               cusum.rec = efp(fit.formula, type="Rec-CUSUM", data=fit.lm$model)
+               plot(cusum.rec, sub=fundName)
+             
+           },
+           "11L"= {
+             ##  CUSUM plot of OLS residuals
+                    
+               cusum.ols = efp(fit.formula, type="OLS-CUSUM", data=fit.lm$model)
+              
+           },
+           "12L"= {
+             ##  CUSUM plot of recursive estimates relative to full sample estimates
+                   
+               cusum.est = efp(fit.formula, type="fluctuation", data=fit.lm$model)
+               plot(cusum.est, functional=NULL, sub=fundName)
+             
+           },
+           "13L"= {
+             ##  rolling regression over 24 month window
+            
+               rollReg <- function(data.z, formula) {
+                 coef(lm(formula, data = as.data.frame(data.z)))  
+               }
+               reg.z = zoo(fit.lm$model, as.Date(rownames(fit.lm$model)))
+               rollReg.z = rollapply(reg.z, FUN=rollReg, fit.formula, width=24, by.column = FALSE, 
+                                     align="right")
+               plot(rollReg.z, main=paste("24-month rolling regression estimates:", fundName, sep=" "))
+            
+           },
+           invisible()
+    )
+    
+    
+    
+  } else {    
+    which.plot<-which.plot[1]
+
+  
+  ##
+  ## 2. Plot selected choices.
+  ##
+ 
+    
+      which.plot <- menu(c("Screeplot of Eigenvalues",
+                          "Factor Returns",
+                            "FM Correlation",
+                            "R square",
+                            "Variance of Residuals",
+                            "Factor Contributions to SD",
+                            "Factor Contributions to ES",
+                            "Factor Contributions to VaR"), title = 
+        "\nMake a plot selection (or 0 to exit):\n")
+      
+      switch(which.plot,
+        "1L" =    {
+  ##
+  ##             1. screeplot.
+  ##
+  if(missing(variables)) {
+    vars <- x$eigen
+    n90 <- which(cumsum(vars)/
+      sum(vars) > 0.9)[1]
+    variables <- 1:max(x$k, min(10, n90))
+  }
+  screeplot(x, variables, cumulative,
+            style, "Screeplot")
+            },
+    "2L" = {
+  ##
+  ##             2. factor returns
+  ##
+  if(missing(variables)) {
+    f.ret <- x$factors
+        }
+    plot.ts(f.ret)
+  
+} ,
+   "3L" = {
+     cov.fm<- factorModelCovariance(t(x$loadings),var(x$factors),x$residVars.vec)    
+     cor.fm = cov2cor(cov.fm)
+     rownames(cor.fm) = colnames(cor.fm)
+     ord <- order(cor.fm[1,])
+     ordered.cor.fm <- cor.fm[ord, ord]
+     plotcorr(ordered.cor.fm, col=cm.colors(11)[5*ordered.cor.fm + 6])  
+   },
+   "4L" ={
+     barplot(x$r2)
+      },
+    "5L" = {
+     barplot(x$residVars.vec)  
+     },
+    "6L" = {
+      cov.factors = var(x$factors)
+      names = colnames(x$asset.ret)
+      factor.sd.decomp.list = list()
+      for (i in names) {
+        factor.sd.decomp.list[[i]] =
+          factorModelSdDecomposition(x$loadings[,i],
+                                     cov.factors, x$residVars.vec[i])
+      }
+      # function to extract contribution to sd from list
+      getCSD = function(x) {
+        x$cr.fm
+      }
+      # extract contributions to SD from list
+      cr.sd = sapply(factor.sd.decomp.list, getCSD)
+      rownames(cr.sd) = c(colnames(x$factors), "residual")
+      # create stacked barchart
+      barplot(cr.sd, main="Factor Contributions to SD",
+              legend.text=T, args.legend=list(x="topleft"),
+              col=c(1:50) )
+    } ,
+    "7L" ={
+      factor.es.decomp.list = list()
+      names = colnames(x$asset.ret)
+      for (i in names) {
+        # check for missing values in fund data
+        idx = which(!is.na(x$asset.ret[,i]))
+        tmpData = cbind(x$asset.ret[idx,i], x$factors,
+                        x$residuals[,i]/sqrt(x$residVars.vec[i]))
+        colnames(tmpData)[c(1,length(tmpData[1,]))] = c(i, "residual")
+        factor.es.decomp.list[[i]] = 
+          factorModelEsDecomposition(tmpData, 
+                                     x$loadings[,i],
+                                     x$residVars.vec[i], tail.prob=0.05)
+      }
+         
+             
+             # stacked bar charts of percent contributions to ES 
+             getCETL = function(x) {
+               x$cES
+             }
+             # report as positive number
+             cr.etl = sapply(factor.es.decomp.list, getCETL)
[TRUNCATED]

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


More information about the Returnanalytics-commits mailing list