[Returnanalytics-commits] r3567 - in pkg/FactorAnalytics: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Dec 4 08:48:23 CET 2014
Author: pragnya
Date: 2014-12-04 08:48:22 +0100 (Thu, 04 Dec 2014)
New Revision: 3567
Added:
pkg/FactorAnalytics/R/plot.sfm.r
pkg/FactorAnalytics/man/plot.sfm.Rd
Modified:
pkg/FactorAnalytics/DESCRIPTION
pkg/FactorAnalytics/NAMESPACE
pkg/FactorAnalytics/R/fitSfm.R
pkg/FactorAnalytics/man/fitSfm.Rd
Log:
Added plot.sfm, corr option for fitSfm
Modified: pkg/FactorAnalytics/DESCRIPTION
===================================================================
--- pkg/FactorAnalytics/DESCRIPTION 2014-12-04 07:46:15 UTC (rev 3566)
+++ pkg/FactorAnalytics/DESCRIPTION 2014-12-04 07:48:22 UTC (rev 3567)
@@ -1,7 +1,7 @@
Package: factorAnalytics
Type: Package
Title: Factor Analytics
-Version: 2.0.4
+Version: 2.0.5
Date: 2014-12-03
Author: Eric Zivot, Sangeetha Srinivasan and Yi-An Chen
Maintainer: Sangeetha Srinivasan <sangee at uw.edu>
Modified: pkg/FactorAnalytics/NAMESPACE
===================================================================
--- pkg/FactorAnalytics/NAMESPACE 2014-12-04 07:46:15 UTC (rev 3566)
+++ pkg/FactorAnalytics/NAMESPACE 2014-12-04 07:48:22 UTC (rev 3567)
@@ -13,6 +13,7 @@
S3method(fmVaRDecomp,sfm)
S3method(fmVaRDecomp,tsfm)
S3method(plot,pafm)
+S3method(plot,sfm)
S3method(plot,tsfm)
S3method(predict,sfm)
S3method(predict,tsfm)
@@ -50,6 +51,7 @@
importFrom(lattice,barchart)
importFrom(lattice,panel.barchart)
importFrom(lattice,panel.grid)
+importFrom(lattice,xyplot)
importFrom(leaps,regsubsets)
importFrom(lmtest,coeftest)
importFrom(lmtest,coeftest.default)
Modified: pkg/FactorAnalytics/R/fitSfm.R
===================================================================
--- pkg/FactorAnalytics/R/fitSfm.R 2014-12-04 07:46:15 UTC (rev 3566)
+++ pkg/FactorAnalytics/R/fitSfm.R 2014-12-04 07:48:22 UTC (rev 3567)
@@ -3,13 +3,15 @@
#' @description Fits a statistical factor model using principal component
#' analysis for one or more asset returns or excess returns. When the number of
#' assets exceeds the number of time periods, APCA (Asymptotic Principal
-#' Component Analysis) is performed. This function is based on the S+FinMetric
-#' function \code{mfactor}. An object of class \code{"sfm"} is returned.
+#' Component Analysis) is performed. An object of class \code{"sfm"} is
+#' returned. This function is based on the S+FinMetric function \code{mfactor}.
#'
#' @details
#' If \code{data} is not of class \code{"xts"}, rownames must provide an
-#' \code{xts} compatible time index. If the data contains missing values,
-#' \code{na.rm} should be set to \code{TRUE} to remove NAs.
+#' \code{xts} compatible time index. Before model fitting, incomplete cases in
+#' \code{data} are removed using \code{\link[stats]{na.omit}}. Also, if
+#' \code{check=TRUE}, a warning is issued if any asset is found to have
+#' identical observations.
#'
#' Let \code{N} be the number of columns or assets and \code{T} be the number
#' of rows or observations. When \code{N < T}, Principal Component Analysis
@@ -27,12 +29,17 @@
#' \code{refine} specifies whether a refinement of the APCA procedure (that may
#' improve efficiency) from Connor and Korajczyk (1988) is to be used.
#'
-#' If \code{check=TRUE}, a warning is issued if any asset is found to have
-#' identical observations.
+#' \code{corr} specifies if the correlation matrix of returns should be used
+#' for finding the principal components instead of the covariance matrix. This
+#' is decided on a case-by-case basis. The variable with the highest variance
+#' dominates the PCA when the covariance matrix is used (but, this may be
+#' justified if a volatile asset is more interesting for some reason and
+#' volatility information shouldn't be discarded). On the other hand, using the
+#' correlation matrix standardizes the variables and makes them comparable,
+#' avoiding penalizing variables with less dispersion.
#'
-#' Note about NAs: Before model fitting, incomplete cases in \code{data} are
-#' removed using \code{\link[stats]{na.omit}}. Otherwise, all observations are
-#' included.
+#' Note: If the median of the 1st principal component is negative, all it's
+#' factor realizations are inverted to enable more meaningful interpretation.
#'
#' @param data vector, matrix, data.frame, xts, timeSeries or zoo object with
#' asset returns. See details.
@@ -46,7 +53,9 @@
#' specified. Default is 0.05.
#' @param check logical; to check if any asset has identical observations.
#' Default is \code{FALSE}.
-#' @param ... arguments passed to other functions.
+#' @param corr logical; whether to use the correlation instead of the covariance
+#' matrix when finding the principal components. Default is \code{FALSE}.
+#' @param ... optional arguments passed to \code{\link[stats]{lm}}.
#'
#' @return fitTsfm returns an object of class \code{"sfm"} for which
#' \code{print}, \code{plot}, \code{predict} and \code{summary} methods exist.
@@ -118,47 +127,45 @@
#' class(sfm.dat)
#' class(sfm.apca.dat)
#'
-#' # pca
+#' # PCA
#' args(fitSfm)
-#' sfm.pca.fit <- fitSfm(sfm.dat, k=2)
-#' class(sfm.pca.fit)
-#' names(sfm.pca.fit)
-#' head(sfm.pca.fit$factors)
-#' head(sfm.pca.fit$loadings)
-#' sfm.pca.fit$r2
-#' sfm.pca.fit$resid.sd
-#' sfm.pca.fit$mimic
+#' fit.pca <- fitSfm(sfm.dat, k=2)
+#' class(fit.pca)
+#' names(fit.pca)
+#' head(fit.pca$factors)
+#' head(fit.pca$loadings)
+#' fit.pca$r2
+#' fit.pca$resid.sd
+#' fit.pca$mimic
#'
-#' # apca with number of factors, k=15
-#' sfm.apca.fit <- fitSfm(sfm.apca.dat, k=15, refine=TRUE)
+#' # APCA with number of factors, k=15
+#' fit.apca <- fitSfm(sfm.apca.dat, k=15, refine=TRUE)
#'
-#' # apca with the Bai & Ng method
-#' sfm.apca.fit.bn <- fitSfm(sfm.apca.dat, k="bn")
+#' # APCA with the Bai & Ng method
+#' fit.apca.bn <- fitSfm(sfm.apca.dat, k="bn")
#'
-#' # apca with the Connor-Korajczyk method
-#' sfm.apca.fit.ck <- fitSfm(sfm.apca.dat, k="ck")
+#' # APCA with the Connor-Korajczyk method
+#' fit.apca.ck <- fitSfm(sfm.apca.dat, k="ck")
#'
#' @importFrom PerformanceAnalytics checkData
#'
#' @export
-fitSfm <- function(data, k=1, max.k=NULL, refine=TRUE, sig=0.05, check=FALSE) {
+fitSfm <- function(data, k=1, max.k=NULL, refine=TRUE, sig=0.05, check=FALSE,
+ corr=FALSE, ...) {
# record the call as an element to be returned
call <- match.call()
- # check input data type and format and coerce to desired type for use
- R.xts <- checkData(data, method="xts")
- R.mat <- coredata(R.xts)
+ # check input data type and coerce to xts; remove NAs
+ R.xts <- na.omit(checkData(data, method="xts"))
- # remove NAs
- R.mat <- na.omit(R.mat)
+ # dim and dimnames of R.mat
+ n <- ncol(R.xts)
+ obs <- nrow(R.xts)
- # dim and dimnames of R.mat
- n <- ncol(R.mat)
- obs <- nrow(R.mat)
- if (is.null(dimnames(data))) {
- dimnames(R.mat) <- list(1:obs, paste("V", 1:n, sep = "."))
+ # assign generic variable names, if they are missing
+ if (is.null(colnames(data))) {
colnames(R.xts) <- paste("V", 1:n, sep = ".")
}
@@ -185,29 +192,32 @@
# check input vailidity or assign default for argument max.k
if (is.null(max.k)) {
max.k <- min(10, obs - 1)
+ } else if (!is.numeric(max.k) || max.k <= 0 || round(max.k) != max.k) {
+ stop("Invalid argument: max.k, the maximum number of factors, must be a
+ positive integer.")
} else if (max.k >= obs) {
stop("Invalid argument: max.k must be less than the number of observations")
}
# check if any asset has identical observations
- temp <- apply(data, 2, range)
- if(any(abs(temp[2, ] - temp[1, ]) < .Machine$single.eps)) {
- warning("Some variables have identical observations.")
+ if (check) {
+ temp <- apply(data, 2, range)
+ if(any(abs(temp[2, ] - temp[1, ]) < .Machine$single.eps)) {
+ warning("Some variables have identical observations.")
+ }
}
# select method to estimate factors based on k and n
# in each case a partial list of return values are obtained
if (n < obs) {
- result <- UsePCA(R.xts=R.xts, R.mat=R.mat, k=k, n=n, obs=obs)
- } else if (k == "ck") {
- result <- UseAPCA_ck(R.xts=R.xts, R.mat=R.mat, max.k=max.k, refine=refine,
- sig=sig, n=n, obs=obs)
- } else if (k == "bn") {
- result <- UseAPCA_bn(R.xts=R.xts, R.mat=R.mat, max.k=max.k, refine=refine,
- n=n, obs=obs)
+ result <- UsePCA(R.xts=R.xts, k=k, corr=corr, ...)
+ } else if (k=="ck") {
+ result <- UseAPCA_ck(R.xts=R.xts, max.k=max.k, refine=refine, sig=sig,
+ corr=corr, ...)
+ } else if (k=="bn") {
+ result <- UseAPCA_bn(R.xts=R.xts, max.k=max.k, refine=refine, corr=corr, ...)
} else {
- result <- UseAPCA(R.xts=R.xts, R.mat=R.mat, k=k, refine=refine, n=n,
- obs=obs)
+ result <- UseAPCA(R.xts=R.xts, k=k, refine=refine, corr=corr, ...)
}
# create list of return values.
@@ -217,14 +227,21 @@
return(result)
}
+
### Principal Component Analysis when N < T
#
-UsePCA <- function(R.xts=R.xts, R.mat=R.mat, k=k, n=n, obs=obs) {
+UsePCA <- function(R.xts=R.xts, k=k, corr=corr, ...) {
+ R.mat <- coredata(R.xts) # TxN
+ n <- ncol(R.mat)
+ obs <- nrow(R.mat)
# demean TxN matrix of returns
R.mat.d <- t(t(R.mat) - colMeans(R.mat))
# NxN return covariance matrix
Omega.N <- crossprod(R.mat.d)/obs
+ if (corr) {
+ Omega.N <- cov2cor(Omega.N)
+ }
# get eigen decomposition
eig.decomp <- eigen(Omega.N, symmetric=TRUE)
eig.val <- eig.decomp$values
@@ -233,10 +250,12 @@
# get TxK factor realizations
f <- R.mat %*% X
colnames(f) <- paste("F", 1:k, sep = ".")
+ # invert 1st principal component if most values are negative
+ if (median(f[,1]) < 0) {f[,1] <- -f[,1]}
# OLS time series regression to get B: NxK matrix of factor loadings
f <- xts(f, index(R.xts))
- asset.fit <- lm(R.xts ~ f)
+ asset.fit <- lm(R.xts ~ f, ...)
B <- t(coef(asset.fit)[-1, , drop=FALSE])
alpha <- coef(asset.fit)[1,]
@@ -252,8 +271,9 @@
mimic <- X / colSums(X)
# assign row and column names
- names(eig.val) = names(r2) = names(resid.sd) = colnames(R.xts)
- colnames(B) <- colnames(f)
+ names(eig.val) <- paste("F", 1:n, sep = ".")
+ names(r2) = names(resid.sd) = colnames(R.xts)
+ colnames(B) = colnames(f)
# return list
list(asset.fit=asset.fit, k=k, factors=f, loadings=B, alpha=alpha, r2=r2,
@@ -264,22 +284,30 @@
### Asymptotic Principal Component Analysis when N >= T
#
-UseAPCA <- function(R.xts=R.xts, R.mat=R.mat, k=k, refine=refine, n=n, obs=obs) {
+UseAPCA <- function(R.xts=R.xts, k=k, refine=refine, corr=corr, ...) {
+ R.mat <- coredata(R.xts) # TxN
+ n <- ncol(R.mat)
+ obs <- nrow(R.mat)
# demean TxN matrix of returns
R.mat.d <- t(t(R.mat) - colMeans(R.mat))
# TxT return covariance matrix
Omega.T <- tcrossprod(R.mat.d)/n
+ if (corr) {
+ Omega.T <- cov2cor(Omega.T)
+ }
# get eigen decomposition
eig.decomp <- eigen(Omega.T, symmetric=TRUE)
eig.val <- eig.decomp$values
# get TxK factor realizations
X <- eig.decomp$vectors[, 1:k, drop=FALSE] # TxK
dimnames(X) <- list(1:obs, paste("F", 1:k, sep = "."))
+ f <- xts(X, index(R.xts))
+ # invert 1st principal component if most values are negative
+ if (median(f[,1]) < 0) {f[,1] <- -f[,1]}
# OLS time series regression to get B: NxK matrix of factor loadings
- f <- xts(X, index(R.xts))
- asset.fit <- lm(R.xts ~ f)
+ asset.fit <- lm(R.xts ~ f, ...)
B <- t(coef(asset.fit)[-1, , drop=FALSE])
alpha <- coef(asset.fit)[1,]
@@ -289,12 +317,16 @@
if (refine) {
R.mat.rescaled <- t(R.mat.d)/resid.sd
Omega.T <- crossprod(R.mat.rescaled)/n
+ if (corr) {
+ Omega.T <- cov2cor(Omega.T)
+ }
eig.decomp <- eigen(Omega.T, symmetric=TRUE)
eig.val <- eig.decomp$values
X <- eig.decomp$vectors[, 1:k, drop=FALSE]
dimnames(X) <- list(1:obs, paste("F", 1:k, sep = "."))
f <- xts(X, index(R.xts))
- asset.fit <- lm(R.xts ~ f)
+ if (median(f[,1]) < 0) {f[,1] <- -f[,1]}
+ asset.fit <- lm(R.xts ~ f, ...)
B <- t(coef(asset.fit)[-1, , drop=FALSE])
alpha <- coef(asset.fit)[1,]
resid.sd <- as.numeric(sapply(X=summary(asset.fit), FUN="[", "sigma"))
@@ -311,9 +343,9 @@
r2 <- as.numeric(sapply(X=summary(asset.fit), FUN="[", "r.squared"))
# assign row and column names
- names(eig.val) = 1:obs
+ names(eig.val) = paste("F", 1:obs, sep = ".")
names(r2) = names(resid.sd) = colnames(R.xts)
- colnames(B) <- colnames(f)
+ colnames(B) = colnames(f)
# return list
list(asset.fit=asset.fit, k=k, factors=f, loadings=B, alpha=alpha, r2=r2,
@@ -324,20 +356,21 @@
### Asymptotic Principal Component Analysis using 'ck' method to determine k
#
-UseAPCA_ck <- function(R.xts=R.xts, R.mat=R.mat, max.k=max.k, refine=refine,
- sig=sig, n=n, obs=obs) {
-
+UseAPCA_ck <- function(R.xts=R.xts, max.k=max.k, refine=refine, sig=sig,
+ corr=corr, ...) {
+ n <- ncol(R.xts)
+ obs <- nrow(R.xts)
idx <- 2*(1:(obs/2))
# dof-adjusted squared residuals for k=1
- fit <- UseAPCA(R.xts=R.xts, R.mat=R.mat, k=1, n=n, obs=obs, refine=refine)
+ fit <- UseAPCA(R.xts=R.xts, k=1, refine=refine, corr=corr, ...)
eps2 <- fit$residuals^2 / (1-2/obs-1/n)
for (k in 2:max.k) {
f <- fit
mu <- rowMeans(eps2[idx-1,,drop=FALSE])
# dof-adjusted squared residuals for k
- fit <- UseAPCA(R.xts=R.xts, R.mat=R.mat, k=k, n=n, obs=obs, refine=refine)
+ fit <- UseAPCA(R.xts=R.xts, k=k, refine=refine, corr=corr, ...)
eps2.star <- fit$residuals^2 / (1-(k+1)/obs-k/n)
mu.star <- rowMeans(eps2.star[idx,,drop=FALSE])
# cross sectional differences in sqd. errors btw odd & even time periods
@@ -352,14 +385,16 @@
### Asymptotic Principal Component Analysis using 'bn' method to determine k
#
-UseAPCA_bn <- function(R.xts=R.xts, R.mat=R.mat, max.k=max.k, refine=refine,
- n=n, obs=obs) {
+UseAPCA_bn <- function(R.xts=R.xts, max.k=max.k, refine=refine, corr=corr, ...) {
+
+ n <- ncol(R.xts)
+ obs <- nrow(R.xts)
# intialize sigma
sigma <- rep(NA, max.k)
for (k in 1:max.k) {
# fit APCA for k factors
- fit <- UseAPCA(R.xts=R.xts, R.mat=R.mat, k=k, n=n, obs=obs, refine=refine)
+ fit <- UseAPCA(R.xts=R.xts, k=k, refine=refine, corr=corr, ...)
# get cross-sectional average of residual variances
sigma[k] <- mean(fit$resid.sd^2)
}
@@ -374,7 +409,7 @@
used.")
}
k <- min(order(PC_p1)[1], order(PC_p2)[1])
- UseAPCA(R.xts=R.xts, R.mat=R.mat, k=k, n=n, obs=obs, refine=refine)
+ UseAPCA(R.xts=R.xts, k=k, refine=refine, corr=corr, ...)
}
Added: pkg/FactorAnalytics/R/plot.sfm.r
===================================================================
--- pkg/FactorAnalytics/R/plot.sfm.r (rev 0)
+++ pkg/FactorAnalytics/R/plot.sfm.r 2014-12-04 07:48:22 UTC (rev 3567)
@@ -0,0 +1,432 @@
+#' @title Plots from a fitted statistical factor model
+#'
+#' @description Generic \code{plot} method for object of class \code{sfm}.
+#' Plots chosen characteristic(s) for one or more assets.
+#'
+#' @details
+#' If the plot type argument is not specified, a menu prompts for user input
+#' and the corresponding plot is output. And, the menu is repeated for
+#' user convenience in plotting multiple characteristics. Selecting '0' from
+#' the menu exits the current \code{plot.sfm} call. Alternately, setting
+#' \code{loop=FALSE} will exit after plotting any one chosen characteristic
+#' without the need for menu selection.
+#'
+#' Group plots are the default. The first \code{n.max} variables and
+#' \code{n.max} factors are plotted depending on the characteristic chosen.
+#'
+#' Individual asset plots are selected by specifying \code{plot.single=TRUE}.
+#' In which case, \code{asset.name} is necessary if multiple assets
+#' were modeled in \code{x}. However, if the \code{fitSfm} object contains only
+#' one asset's factor model fit, \code{plot.sfm} can infer this automatically,
+#' without user input.
+#'
+#' @param x an object of class \code{sfm} produced by \code{fitSfm}.
+#' @param which.plot.group a number to indicate the type of group plot for
+#' multiple assets. If \code{NULL} (default), the following menu appears: \cr
+#' 1 = Screeplot of eigenvalues, \cr
+#' 2 = Time series plot of estimated factors, \cr
+#' 3 = Estimated factor loadings, \cr
+#' 4 = Histogram of R-squared, \cr
+#' 5 = Histogram of Residual volatility,\cr
+#' 6 = Factor model residual correlation \cr
+#' 7 = Factor model correlation,\cr
+#' 8 = Factor contribution to SD,\cr
+#' 9 = Factor contribution to ES,\cr
+#' 10 = Factor contribution to VaR
+#' @param k.max maximum number of factors to add per plot device for group
+#' plots. Default is 6.
+#' @param n.max maximum number of variables to add per plot device for group
+#' plots. Default is 10.
+#' @param plot.single logical; If \code{TRUE} plots the characteristics of an
+#' individual asset's factor model. The type of plot is given by
+#' \code{which.plot.single}. Default is \code{FALSE}.
+#' @param asset.name name of the individual asset to be plotted. Is necessary
+#' if \code{x} contains multiple asset fits and \code{plot.single=TRUE}.
+#' @param which.plot.single a number to indicate the type of group plot for an
+#' individual asset. If \code{NULL} (default), the following menu appears: \cr
+#' 1 = Time series plot of actual and fitted asset returns,\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 test-Recursive residuals,\cr
+#' 11 = CUSUM test-OLS residuals,\cr
+#' 12 = Recursive estimates (RE) test of OLS regression coefficients,\cr
+#' 13 = Rolling estimates over a 24-period observation window
+#' @param colorset color palette to use for all the plots. Default is
+#' \code{c(1:12)}. The 1st element will be used for individual time series
+#' plots or the 1st series plotted, the 2nd element for the 2nd object in the
+#' plot and so on.
+#' @param legend.loc places a legend into one of nine locations on the chart:
+#' "bottomright", "bottom", "bottomleft", "left", "topleft", "top", "topright",
+#' "right", or "center". Default is "bottomright". Use \code{legend.loc=NULL}
+#' to suppress the legend.
+#' @param las one of {0, 1, 2, 3} to set the direction of axis labels, same as
+#' in \code{plot}. Default here is 1.
+#' @param VaR.method a method for computing VaR; one of "modified", "gaussian",
+#' "historical" or "kernel". VaR is computed using
+#' \code{\link[PerformanceAnalytics]{VaR}}. Default is "historical".
+#' @param eig.max scalar in (0,1] for limiting the screeplot to factors that
+#' explain a given percent of the variance. Default is 0.9.
+#' @param cum.var logical; If \code{TRUE}, the cumulative fraction of the
+#' variance is printed above each bar in the screeplot of eigenvalues. Default
+#' is \code{TRUE}.
+#' @param loop logical to indicate if the plot menu should be repeated. Default
+#' is \code{TRUE}.
+#' @param ... further arguments to be passed to other plotting functions.
+#'
+#' @author Eric Zivot, Sangeetha Srinivasan and Yi-An Chen
+#'
+#' @seealso \code{\link{fitSfm}} and \code{\link{summary.sfm}} for details
+#' about the time series factor model fit, extractor functions and summary
+#' statistics.
+#'
+#' \code{\link[strucchange]{efp}} for CUSUM tests.
+#'
+#' \code{\link[xts]{plot.xts}},
+#' \code{\link[PerformanceAnalytics]{chart.TimeSeries}},
+#' \code{\link[PerformanceAnalytics]{chart.ACFplus}},
+#' \code{\link[PerformanceAnalytics]{chart.Histogram}},
+#' \code{\link[PerformanceAnalytics]{chart.QQPlot}},
+#' \code{\link[graphics]{barplot}}, \code{\link[lattice]{barchart}} and
+#' \code{\link[corrplot]{corrplot}} for plotting methods used.
+#'
+#' \code{\link{fmSdDecomp}}, \code{\link{fmEsDecomp}},
+#' \code{\link{fmVaRDecomp}} for factor model risk measures.
+#'
+#' @examples
+#'
+#' # load data from the database
+#' data(stat.fm.data)
+#'
+#' # APCA with number of factors, k=15
+#' fit.apca <- fitSfm(sfm.apca.dat, k=15, refine=TRUE)
+#'
+#' # group plot(default); select type from menu prompt
+#' # menu is auto-looped to get multiple types of plots based on the same fit
+#' # plot(fit.apca)
+#'
+#' # plot the factor betas of 1st 4 assets fitted above
+#' # loop disabled to get one type of plot without interative menu
+#' plot(fit.apca, n.max=4, which.plot.group=3, loop=FALSE)
+#'
+#' # plot factor model return correlation; angular order of the eigenvectors
+#' plot(fit.apca, which.plot.group=7, loop=FALSE,
+#' order="AOE", method="ellipse", tl.pos = "d")
+#'
+#' # histogram of residuals from an individual asset's factor model fit
+#' plot(fit.apca, plot.single=TRUE, asset.name="AFL", which.plot.single=8,
+#' loop=FALSE)
+#'
+#' @importFrom PerformanceAnalytics chart.TimeSeries chart.ACFplus
+#' chart.Histogram chart.QQPlot
+#' @importFrom lattice barchart xyplot panel.barchart panel.grid
+#' @importFrom corrplot corrplot
+#' @importFrom strucchange efp
+#'
+#' @method plot sfm
+#' @export
+
+plot.sfm <- function(x, which.plot.group=NULL, k.max=6, n.max=10,
+ plot.single=FALSE, asset.name, which.plot.single=NULL,
+ colorset=(1:12), legend.loc="topleft", las=1,
+ VaR.method="historical", cum.var=TRUE, eig.max=0.9,
+ loop=TRUE, ...) {
+
+ if (plot.single==TRUE) {
+
+ if (missing(asset.name) && length(x$asset.names)>1) {
+ stop("Missing input: 'asset.name' is required if plot.single is TRUE and
+ the factor model fits multiple assets.")
+ } else if (length(x$asset.names)==1) {
+ i <- x$asset.names[1]
+ } else {
+ i <- asset.name
+ }
+ # extract info from the fitSfm object
+ plotData <- merge.xts(x$data[,i], fitted(x)[,i])
+ colnames(plotData) <- c("Actual","Fitted")
+ Residuals <- residuals(x)[,i]
+ fit <- lm(x$data[,i] ~ x$factors)
+
+ # plot selection
+ repeat {
+ if (is.null(which.plot.single)) {
+ which.plot.single <-
+ menu(c("Time series plot of actual and fitted asset returns",
+ "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 test-Recursive residuals",
+ "CUSUM test-OLS residuals",
+ "Recursive estimates (RE) test of OLS regression coefficients",
+ "Rolling estimates over a 24-period observation window"),
+ title="\nMake a plot selection (or 0 to exit):")
+ }
+
+ par(las=las) # default horizontal axis labels
+
+ switch(which.plot.single,
+ "1L" = {
+ ## time series plot of actual and fitted asset returns
+ chart.TimeSeries(plotData, main=paste("Returns:",i),
+ colorset=colorset, xlab="",
+ ylab="Actual and fitted asset returns",
+ legend.loc=legend.loc, pch=NULL, las=las, ...)
+ }, "2L" = {
+ ## time series plot of residuals with standard error bands
+ chart.TimeSeries(Residuals, main=paste("Residuals:",i),
+ colorset=colorset, xlab="", ylab="Residuals",
+ lwd=2, lty="solid", las=las, ...)
+ abline(h=1.96*x$resid.sd[i], lwd=2, lty="dotted", col="red")
+ abline(h=-1.96*x$resid.sd[i], lwd=2, lty="dotted", col="red")
+ legend(x=legend.loc, lty=c("solid","dotted"),
+ col=c(colorset[1],"red"), lwd=2,
+ legend=c("Residuals",expression("\u00b1 1.96"*sigma)))
+ }, "3L" = {
+ ## time series plot of squared residuals
+ chart.TimeSeries(Residuals^2, colorset=colorset, xlab="",
+ ylab=" Squared Residuals",
+ main=paste("Squared Residuals:",i),
+ legend.loc=legend.loc, pch=NULL, las=las, ...)
+ }, "4L" = {
+ ## time series plot of absolute residuals
+ chart.TimeSeries(abs(Residuals), colorset=colorset, xlab="",
+ ylab="Absolute Residuals",
+ main=paste("Absolute Residuals:",i),
+ legend.loc=legend.loc, pch=NULL, las=las, ...)
+ }, "5L" = {
+ ## SACF and PACF of residuals
+ chart.ACFplus(Residuals, col=colorset[1],
+ main=paste("SACF & PACF - Residuals:",i), ...)
+ }, "6L" = {
+ ## SACF and PACF of squared residuals
+ chart.ACFplus(Residuals^2, col=colorset[1], ...,
+ main=paste("SACF & PACF - Squared residuals:",i))
+ }, "7L" = {
+ ## SACF and PACF of absolute residuals
+ chart.ACFplus(abs(Residuals), col=colorset[1], ...,
+ main=paste("SACF & PACF - Absolute Residuals:",i))
+ }, "8L" = {
+ ## histogram of residuals with normal curve overlayed
+ methods <- c("add.density","add.normal","add.rug")
+ chart.Histogram(Residuals, xlab="Return residuals",
+ methods=methods, colorset=colorset,
+ main=paste("Histogram of Residuals:",i), ...)
+ }, "9L" = {
+ ## normal qq-plot of residuals
+ chart.QQPlot(Residuals, envelope=0.95, col=colorset,
+ main=paste("QQ-plot of Residuals:",i), ...)
+ legend(x=legend.loc, col="red", lty="dotted", lwd=1,
+ legend=c("0.95 confidence envelope"))
+ }, "10L" = {
+ ## Recursive CUSUM test
+ cusum.rec <- efp(formula(fit), type="Rec-CUSUM", data=fit$model)
+ plot(cusum.rec, main=paste("Recursive CUSUM test:",i), las=las,
+ col=colorset, ...)
+ }, "11L" = {
+ ## OLS-based CUSUM test
+ cusum.ols <- efp(formula(fit), type="OLS-CUSUM", data=fit$model)
+ plot(cusum.ols, main=paste("OLS-based CUSUM test:",i), las=las,
+ col=colorset, ...)
+ }, "12L" = {
+ ## Recursive estimates (RE) test of OLS regression coefficients
+ cusum.est <- efp(formula(fit), type="RE", data=fit$model)
+ plot(cusum.est, functional=NULL, col=colorset, las=0,
+ main=paste("RE test (Recursive estimates test):",i), ...)
+ }, "13L" = {
+ ## Rolling estimates over 24-period observation window
+ rollReg <- function(data.z, formula) {
+ coef(lm(formula, data=as.data.frame(data.z)))
+ }
+ reg.z <- zoo(fit$model, as.Date(rownames(fit$model)))
+ rollReg.z <- rollapply(reg.z, FUN=rollReg, formula(fit),
+ width=24, by.column=FALSE, align="right")
+ par(las=0)
+ plot(rollReg.z, ..., las=las,
+ main=paste("Rolling estimates (24-period obs window):",i))
+ par(las=las)
+ },
+ invisible()
+ )
+ # repeat menu if user didn't choose to exit from the plot options
+ if (which.plot.single==0 || loop==FALSE) {break}
+ else {which.plot.single=NULL}
+ }
+ } else { # start of group asset plots
+
+ # plot selection
+ repeat {
+ if (is.null(which.plot.group)) {
+ which.plot.group <-
+ menu(c("Screeplot of eigenvalues",
+ "Time series plot of estimated factors",
+ "Estimated factor loadings",
+ "Histogram of R-squared",
+ "Histogram of Residual volatility",
+ "Factor model residual correlation",
+ "Factor model return correlation",
+ "Factor contribution to SD",
+ "Factor contribution to ES",
+ "Factor contribution to VaR"),
+ title="\nMake a plot selection (or 0 to exit):")
+ }
+
+ par(las=las) # default horizontal axis labels
+ k <- x$k
+ n <- nrow(x$loadings)
+
+ switch(which.plot.group,
+ "1L" = {
+ ## Screeplot of eigenvalues
+ cumv <- cumsum(x$eigen)/sum(x$eigen)
+ limit <- length(cumv[cumv<eig.max]) + 1
+ eig.pct <- (x$eigen/sum(x$eigen))[1:limit]
+ scree <- barplot(eig.pct, main="Screeplot of eigenvalues",
+ ylab="Proportion of Variance", col="darkblue",
+ ylim=c(0, 1.1*max(eig.pct)), las=las, ...)
+ if (cum.var) {
+ text(scree, eig.pct, label=round(cumv,3), pos=3, cex=0.75)
+ }
+ },
+ "2L" = {
+ ## Time series plot of estimated factors
+ if (k > k.max) {
+ cat(paste("Displaying only the first", k.max,"factors, as the number of factors > 'k.max' =", k.max))
+ k <- k.max
+ }
+ plot(
+ xyplot(x$factors[,1:k],type=c("l","g"),xlab="",
+ scales=list(y=list(rot=0)), strip.left=TRUE, strip=FALSE)
+ )
+ },
+ "3L" = {
+ ## Estimated factor loadings
+ if (k > k.max) {
+ cat(paste("Displaying only the first", k.max,"factors, as the number of factors > 'k.max' =", k.max))
+ k <- k.max
+ }
+ if (n > n.max) {
+ cat(paste("Displaying only the first", n.max,"variables, as the number of variables > 'n.max' =", n.max))
+ n <- n.max
+ }
+ par(mfrow=c(ceiling(k/2),2))
+ for (i in 1:k) {
+ main=paste("Beta values for ", colnames(x$loadings)[i])
+ barplot(x$loadings[1:n,i], main=main, names.arg=x$asset.names[1:n],
+ col="darkblue", las=las, horiz=TRUE, ...)
+ abline(v=0, lwd=1, lty=1, col=1)
+ }
+ par(mfrow=c(1,1))
+ },
+ "4L" ={
+ ## Histogram of R-squared
+ methods <- c("add.density","add.rug")
+ chart.Histogram(x$r2, xlab="R-squared",
+ methods=methods, colorset=colorset,
+ main=paste("Histogram of R-squared values"), ...)
+ },
+ "5L" = {
+ ## Histogram of Residual Volatility
+ methods <- c("add.density","add.rug")
+ chart.Histogram(x$resid.sd, xlab="Residual volatility",
+ methods=methods, colorset=colorset,
+ main=paste("Histogram of Residual volatilities"), ...)
+ },
+ "6L" = {
+ ## Factor Model Residual Correlation
+ if (n > n.max) {
+ cat(paste("Displaying only the first", n.max,"variables, as the number of variables > 'n.max' =", n.max))
+ n <- n.max
+ }
+ cor.resid <- cor(residuals(x)[,1:n], use="pairwise.complete.obs")
+ corrplot(cor.resid, ...)
+ # mtext("pairwise complete obs", line=0.5)
+ },
+ "7L" = {
+ ## Factor Model Return Correlation
+ if (n > n.max) {
+ cat(paste("Displaying only the first", n.max,"variables, as the number of variables > 'n.max' =", n.max))
+ n <- n.max
+ }
+ cor.fm <- cov2cor(fmCov(x))[1:n,1:n]
+ corrplot(cor.fm, ...)
+ # mtext("pairwise complete obs", line=0.5)
+ },
+ "8L" = {
+ ## Factor Percentage Contribution to SD
+ if (k > k.max) {
+ cat(paste("Displaying only the first", k.max,"factors, as the number of factors > 'k.max' =", k.max))
+ k <- k.max
+ }
+ if (n > n.max) {
+ cat(paste("Displaying only the first", n.max,"variables, as the number of variables > 'n.max' =", n.max))
+ n <- n.max
+ }
+ pcSd.fm <- fmSdDecomp(x)$pcSd[1:n,1:k]
+ plot(
+ barchart(pcSd.fm, main="Factor % Contribution to SD", xlab="",
+ auto.key=list(space="bottom",columns=3,
+ points=FALSE,rectangles=TRUE),
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/returnanalytics -r 3567
More information about the Returnanalytics-commits
mailing list