#' plot StatFactorModel object.
#'
#' Generic function of plot method for fitStatisticFactorModel. Either plot all
#' assets or choose a single asset to plot.
#'
#'
#' @param x fit object created by \code{fitStatisticalFactorModel}.
#' @param variables Optional. An integer vector shows which variables are to
#' be plotted. The default is to plot all the variables, or the number of
#' variables explaining 90 percent of the variance, whatever is bigger.
#' @param cumulative Logical flag: if \code{TRUE}, the cumulative fraction of the
#' variance is printed above each bar in the plot.
#' @param style Charater. bar or lines can be chosen.
#' @param which.plot integer indicating which plot to create: "none" will
#' create a menu to choose. Defualt is none.\cr
#' 1 = "Screeplot of Eigenvalues", \cr
#' 2 = "Factor returns", \cr
#' 3 = "FM Correlation", \cr
#' 4 = "R square",\cr
#' 5 = "Variance of Residuals", \cr
#' 6 = "Factor Contributions to SD", \cr
#' 7 = "Factor Contributions to ES", \cr
#' 8 = "Factor Contributions to VaR" \cr
#' @param hgrid Logic flag. Whether to plot horizontal grid or not. Defualt is
#' \code{FALSE}.
#' @param vgrid Logic flag. Whether to plot vertical grid or not. Defualt is \code{FALSE}.
#' @param plot.single Plot a single asset of lm class. Defualt is \code{FALSE}.
#' @param asset.name Name of the asset to be plotted.
#' @param which.plot.single integer indicating which plot to create: "none"
#' will create a menu to choose. Defualt is none. \cr
#' 1 = time series plot of actual and fitted values, \cr
#' 2 = time series plot of residuals with standard error bands,\cr
#' 3 = time series plot of squared residuals, \cr
#' 4 = time series plot of absolute residuals, \cr
#' 5 = SACF and PACF of residuals, \cr
#' 6 = SACF and PACF of squared residuals, \cr
#' 7 = SACF and PACF of absolute residuals, \cr
#' 8 = histogram of residuals with normal curve overlayed, \cr
#' 9 = normal qq-plot of residuals,\cr
#' 10= CUSUM plot of recursive residuals,\cr
#' 11= CUSUM plot of OLS residuals,\cr
#' 12= CUSUM plot of recursive estimates relative to full sample estimates,\cr
#' 13= rolling estimates over 24 month window.
#' @param max.show Maximum assets to plot. Default is 6.
#' @param VaR.method Character, method for computing VaR. Valid choices are
#' either "modified","gaussian","historical", "kernel". computation is done with the \code{VaR}
#' in the PerformanceAnalytics package. Default is "historical".
#' @param ... Other variables for barplot method.
#' @author Eric Zivot and Yi-An Chen.
#' @examples
#'
#' \dontrun{
#' # load data for fitStatisticalFactorModel.r
#' # data from finmetric berndt.dat and folio.dat
#'
#' data(stat.fm.data)
#' # pca
#' sfm.pca.fit <- fitStatisticalFactorModel(sfm.dat,k=10)
#' args(plot.StatFactorModel)
#' # plot all
#' plot(sfm.pca.fit)
#' # plot single asset
#' plot(sfm.pca.fit,plot.single=TRUE,asset.name="CITCRP")
#' }
#' @method plot StatFactorModel
#' @export
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, asset.name,
which.plot.single=c("none","1L","2L","3L","4L","5L","6L",
"7L","8L","9L","10L","11L","12L","13L"),
max.show=6, VaR.method = "historical",...)
{
require(strucchange)
require(ellipse)
#
# beginning of funciton screenplot
#
screeplot<-
function(mf, variables, cumulative = TRUE, style = "bar", main = "", ...)
{
vars <- mf$eigen
if(missing(variables)) {
variables <- 1:mf$k
}
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]/sum(vars),
names = paste("F", variables,sep = "."),
main = main, ylab = "Percentage of Variances", ...)
}
else {
loc <- 1:length(variables)
plot(loc, vars[variables]/sum(vars), type = "b", axes = F, main = main,
ylab = "Percentage of 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.
## asset.name charater. The name of the single asset to be ploted.
## 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]
if (which.plot.single=="none")
# pca method
if ( dim(x$asset.ret)[1] > dim(x$asset.ret)[2] ) {
fit.lm = x$asset.fit[[asset.name]]
## exact information from lm object
factorNames = colnames(fit.lm$model)[-1]
fit.formula = as.formula(paste(asset.name,"~", paste(factorNames, collapse="+"), sep=" "))
#Date = try(as.Date(names(residuals(fit.lm))))
#Date = try(as.yearmon(names(residuals(fit.lm)),"%b %Y"))
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)
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=asset.name, ylab="Monthly performance", lwd=2, col="black")
lines(fitted.z, lwd=2, col="red")
abline(h=0)
legend(x="bottomleft", legend=c("Actual", "Fitted"), lwd=2, col=c("black","red"))
},
"2L" = {
## time series plot of residuals with standard error bands
plot(residuals.z, main=asset.name, 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=asset.name, 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=asset.name, 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: ", asset.name, sep=""))
},
"6L" = {
## SACF and PACF of squared residuals
chart.ACFplus(residuals.z^2, main=paste("Residuals^2: ", asset.name, sep=""))
},
"7L" = {
## SACF and PACF of absolute residuals
chart.ACFplus(abs(residuals.z), main=paste("|Residuals|: ", asset.name, sep=""))
},
"8L" = {
## histogram of residuals with normal curve overlayed
chart.Histogram(residuals.z, methods="add.normal", main=paste("Residuals: ", asset.name, sep=""))
},
"9L" = {
## normal qq-plot of residuals
chart.QQPlot(residuals.z, envelope=0.95, main=paste("Residuals: ", asset.name, sep=""))
},
"10L"= {
## CUSUM plot of recursive residuals
cusum.rec = efp(fit.formula, type="Rec-CUSUM", data=fit.lm$model)
plot(cusum.rec, sub=asset.name)
},
"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=asset.name)
},
"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:", asset.name, sep=" "))
},
invisible()
)
} else { #apca method
dates <- names(x$data[,asset.name])
actual.z <- zoo(x$asset.ret[,asset.name],as.Date(dates))
residuals.z <- zoo(x$residuals,as.Date(dates))
fitted.z <- actual.z - residuals.z
t <- length(dates)
k <- x$k
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"),
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[,asset.name], main=asset.name, ylab="Monthly performance", lwd=2, col="black")
lines(fitted.z[,asset.name], lwd=2, col="red")
abline(h=0)
legend(x="bottomleft", legend=c("Actual", "Fitted"), lwd=2, col=c("black","red"))
},
"2L"={
# "time series plot of residuals with standard error bands"
plot(residuals.z[,asset.name], main=asset.name, ylab="Monthly performance", lwd=2, col="black")
abline(h=0)
sigma = (sum(residuals.z[,asset.name]^2)*(t-k)^-1)^(1/2)
abline(h=2*sigma, lwd=2, lty="dotted", col="red")
abline(h=-2*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[,asset.name]^2, main=asset.name, 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[,asset.name]), main=asset.name, 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[,asset.name], main=paste("Residuals: ", asset.name, sep=""))
},
"6L" = {
## SACF and PACF of squared residuals
chart.ACFplus(residuals.z[,asset.name]^2, main=paste("Residuals^2: ", asset.name, sep=""))
},
"7L" = {
## SACF and PACF of absolute residuals
chart.ACFplus(abs(residuals.z[,asset.name]), main=paste("|Residuals|: ", asset.name, sep=""))
},
"8L" = {
## histogram of residuals with normal curve overlayed
chart.Histogram(residuals.z[,asset.name], methods="add.normal", main=paste("Residuals: ", asset.name, sep=""))
},
"9L" = {
## normal qq-plot of residuals
chart.QQPlot(residuals.z[,asset.name], envelope=0.95, main=paste("Residuals: ", asset.name, sep=""))
},
invisible() )
}
} else {
which.plot<-which.plot[1]
##
## 2. Plot selected choices.
##
if(which.plot=='none')
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
variables <- 1:x$k
}
screeplot(x, variables, cumulative,
style, "Screeplot of Eigenvalues")
},
"2L" = {
##
## 2. factor returns
##
if(missing(variables)) {
f.ret <- x$factors
}
plot.zoo(f.ret)
} ,
"3L" = {
cov.fm<- factorModelCovariance(t(x$loadings),var(x$factors),
x$resid.variance)
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[(1:max.show),(1:max.show)], col=cm.colors(11)[5*ordered.cor.fm + 6])
},
"4L" ={
barplot(x$r2[1:max.show])
},
"5L" = {
barplot(x$resid.variance[1:max.show])
},
"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$resid.variance[i])
}
# function to extract component contribution to sd from list
getCSD = function(x) {
x$cSd.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[,(1:max.show)], main="Factor Contributions to SD",...)
# legend.text=T, args.legend=list(x="topright"))
} ,
"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$resid.variance[i]))
colnames(tmpData)[c(1,length(tmpData[1,]))] = c(i, "residual")
factor.es.decomp.list[[i]] =
factorModelEsDecomposition(tmpData,
x$loadings[,i],
x$resid.variance[i], tail.prob=0.05,VaR.method=VaR.method)
}
# stacked bar charts of component contributions to ES
getCETL = function(x) {
x$cES
}
# report as positive number
cr.etl = sapply(factor.es.decomp.list, getCETL)
rownames(cr.etl) = c(colnames(x$factors), "residual")
barplot(cr.etl[,(1:max.show)], main="Factor Contributions to ES",...)
# legend.text=T, args.legend=list(x="topright") )
},
"8L" = {
factor.VaR.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$resid.variance[i]))
colnames(tmpData)[c(1,length(tmpData[1,]))] = c(i, "residual")
factor.VaR.decomp.list[[i]] =
factorModelVaRDecomposition(tmpData,
x$loadings[,i],
x$resid.variance[i], tail.prob=0.05,VaR.method=VaR.method)
}
# stacked bar charts of component contributions to VaR
getCVaR = function(x) {
x$cVaR.fm
}
# report as positive number
cr.var = sapply(factor.VaR.decomp.list, getCVaR)
rownames(cr.var) = c(colnames(x$factors), "residual")
barplot(cr.var[,(1:max.show)], main="Factor Contributions to VaR",...)
# legend.text=T, args.legend=list(x="topright"))
}, invisible()
)
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.