[Uwgarp-commits] r167 - in pkg/GARPFRM: demo sandbox vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Apr 16 01:12:17 CEST 2014
Author: rossbennett34
Date: 2014-04-16 01:12:16 +0200 (Wed, 16 Apr 2014)
New Revision: 167
Added:
pkg/GARPFRM/demo/CAPM_TF.R
pkg/GARPFRM/demo/DataAccess.R
pkg/GARPFRM/demo/DelineatingEfficientPortfolios.R
pkg/GARPFRM/demo/EstimatingVolatilitiesCorrelation.R
pkg/GARPFRM/demo/MonteCarloMethods.R
pkg/GARPFRM/demo/PerformanceMeasures.R
pkg/GARPFRM/demo/QuantifyingVolatilityVaRModels.R
pkg/GARPFRM/demo/QuantitativeAnalysisBasics.R
pkg/GARPFRM/demo/WebApplications.R
pkg/GARPFRM/sandbox/vignette_to_demo.R
Removed:
pkg/GARPFRM/vignettes/PA_vignette.pdf
Log:
Adding script to pull R chunks from vignettes and write to demo files. Adding demos based on R code from vignettes
Added: pkg/GARPFRM/demo/CAPM_TF.R
===================================================================
--- pkg/GARPFRM/demo/CAPM_TF.R (rev 0)
+++ pkg/GARPFRM/demo/CAPM_TF.R 2014-04-15 23:12:16 UTC (rev 167)
@@ -0,0 +1,82 @@
+
+
+# 'Load the GARPFRM package and CRSP dataset for CAPM analysis.
+suppressMessages(library(GARPFRM))
+options(digits=3)
+data(crsp.short)
+
+stock.df <- largecap.ts[, 1:20]
+mrkt <- largecap.ts[, "market"]
+rfr <- largecap.ts[, "t90"]
+
+# Plot first four stocks from
+plot.zoo(stock.df[,1:4], main="First Four Large Cap Returns")
+
+
+
+# Illustrate the type of data being analzyed: start-end dates.
+start(stock.df[,1:4])
+end(stock.df[,1:4])
+# Count the number of rows: sample size.
+nrow(stock.df)
+
+
+
+# Excess Returns initialized before utilizing in CAPM
+exReturns <- Return.excess(stock.df, rfr)
+colnames(exReturns)= c(colnames(stock.df))
+
+
+
+# Univariate CAPM
+uv <- CAPM(exReturns[,1], mrkt)
+getStatistics(uv)
+
+# Plot data with regression line
+plot(uv)
+
+
+
+# MLM CAPM for AMAT, AMGN, and CAT
+mlm <- CAPM(exReturns[,1:3], mrkt)
+getStatistics(mlm)
+
+# Plot data with regression line
+plot(mlm)
+
+
+
+# For uv
+getBetas(uv)
+getAlphas(uv)
+hypTest(uv, significanceLevel=0.05)
+# For mlm
+getBetas(mlm)
+getAlphas(mlm)
+hypTest(mlm, significanceLevel=0.05)
+
+
+
+# MLM CAPM
+mlm <- CAPM(exReturns[,], mrkt)
+
+# Plot expected returns versus betas
+chartSML(mlm)
+
+
+
+# Load FED consumption data: CONS
+data(consumption)
+
+# Convert to yearmon index and align consumption and mrkt
+consumption <- xts(consumption, as.yearmon(index(consumption)))
+mrkt <- xts(mrkt, as.yearmon(index(mrkt)))
+consumption <- consumption[index(mrkt)]
+
+capm.cons = CAPM(consumption, mrkt)
+coef(summary(capm.cons))
+
+# Plot data with regression line
+plot(capm.cons)
+
+
Added: pkg/GARPFRM/demo/DataAccess.R
===================================================================
--- pkg/GARPFRM/demo/DataAccess.R (rev 0)
+++ pkg/GARPFRM/demo/DataAccess.R 2014-04-15 23:12:16 UTC (rev 167)
@@ -0,0 +1,22 @@
+
+
+library('knitr')
+opts_chunk$set(message=FALSE, fig.path='figures/', fig.align='center', fig.width=4, fig.height=3, fig.keep='last', dev.args=list(pointsize=8))
+options(width=80)
+
+
+
+library(quantmod)
+args(getSymbols)
+
+
+
+getSymbols('^GSPC')
+chart_Series(GSPC)
+
+
+
+getSymbols('DGS3MO',src='FRED')
+plot(DGS3MO,main="3-Month Treasury Constant Maturity Rate",cex.main=0.75)
+
+
Added: pkg/GARPFRM/demo/DelineatingEfficientPortfolios.R
===================================================================
--- pkg/GARPFRM/demo/DelineatingEfficientPortfolios.R (rev 0)
+++ pkg/GARPFRM/demo/DelineatingEfficientPortfolios.R 2014-04-15 23:12:16 UTC (rev 167)
@@ -0,0 +1,291 @@
+
+
+library(knitr)
+opts_chunk$set(tidy=FALSE, warning=FALSE, fig.width=5, fig.height=5)
+
+
+
+suppressPackageStartupMessages(library(GARPFRM))
+
+# Colonel Motors expected return
+R_C <- 0.14
+
+# Colonel Motors standard deviation
+sigma_C <- 0.06
+
+# Separated Edison expected return
+R_S <- 0.08
+
+# Separated Edison standard deviation
+sigma_S <- 0.03
+
+# Correlation coefficient
+rho <- 1
+
+# Create a vector of portfolio weights
+X_C <- seq(from=0, to=1, by=0.2)
+
+# Calculate the portfolio expected return
+R_P <- portReturnTwoAsset(R_C, R_S, X_C)
+
+# Calculate the portfolio standard deviation
+sigma_P <- portSDTwoAsset(R_C, R_S, X_C, sigma_C, sigma_S, rho)
+
+# Combine the portfolio returns and standard deviations in a data.frame object.
+# Note that this is the transpose t() of the data.frame object to match the
+# layout of the tables in the FRM book.
+df <- t(data.frame(R_P=R_P, sigma_P=sigma_P))
+colnames(df) <- X_C
+df
+
+
+
+# Create and plot the efficient frontier object
+ef <- efficientFrontierTwoAsset(R_C, R_S, sigma_C, sigma_S, rho)
+plot(ef)
+
+
+
+# Correlation coefficient
+rho <- -1
+
+# Calculate the portfolio expected return
+R_P <- portReturnTwoAsset(R_C, R_S, X_C)
+
+# Calculate the portfolio standard deviation
+sigma_P <- portSDTwoAsset(R_C, R_S, X_C, sigma_C, sigma_S, rho)
+
+# Combine the portfolio returns and standard deviations in a data.frame object.
+# Note that this is the transpose t() of the data.frame object to match the
+# layout of the tables in the FRM book.
+df <- t(data.frame(R_P=R_P, sigma_P=sigma_P))
+colnames(df) <- X_C
+df
+
+
+
+# Create and plot the efficient frontier object
+ef <- efficientFrontierTwoAsset(R_C, R_S, sigma_C, sigma_S, rho)
+plot(ef)
+
+
+
+# Correlation coefficient
+rho <- 0
+
+# Calculate the portfolio expected return
+R_P <- portReturnTwoAsset(R_C, R_S, X_C)
+
+# Calculate the portfolio standard deviation
+sigma_P <- portSDTwoAsset(R_C, R_S, X_C, sigma_C, sigma_S, rho)
+
+# Combine the portfolio returns and standard deviations in a data.frame object.
+# Note that this is the transpose t() of the data.frame object to match the
+# layout of the tables in the FRM book.
+df <- t(data.frame(R_P=R_P, sigma_P=sigma_P))
+colnames(df) <- X_C
+df
+
+
+
+# Create and plot the efficient frontier object
+ef <- efficientFrontierTwoAsset(R_C, R_S, sigma_C, sigma_S, rho)
+plot(ef)
+
+
+
+
+minRisk <- function(R_A, R_B, sigma_A, sigma_B, rho){
+ top_term <- sigma_B^2 - sigma_A * sigma_B * rho
+ bottom_term <- sigma_A^2 + sigma_B^2 - 2 * sigma_A * sigma_B * rho
+ # x is the fraction of asset A
+ x <- top_term / bottom_term
+ # calculate the weights for the minimum risk portfolio
+ weights <- c(x, 1-x)
+ names(weights) <- c("Asset_A", "Asset_B")
+ # calculate the portfolio return and standard deviation of the minimum
+ # risk portfolio
+ port_ret <- portReturnTwoAsset(R_A, R_B, x)
+ port_sd <- portSDTwoAsset(R_A, R_B, x, sigma_A, sigma_B, rho)
+ return(list(weights=weights,
+ portfolio_return=port_ret,
+ portfolio_sd=port_sd))
+}
+
+
+
+minRisk(R_C, R_S, sigma_C, sigma_S, rho)
+
+
+
+# Correlation coefficient
+rho <- 0.5
+
+# Calculate the portfolio expected return
+R_P <- portReturnTwoAsset(R_C, R_S, X_C)
+
+# Calculate the portfolio standard deviation
+sigma_P <- portSDTwoAsset(R_C, R_S, X_C, sigma_C, sigma_S, rho)
+
+# Combine the portfolio returns and standard deviations in a data.frame object.
+# Note that this is the transpose t() of the data.frame object to match the
+# layout of the tables in the FRM book.
+df <- t(data.frame(R_P=R_P, sigma_P=sigma_P))
+colnames(df) <- X_C
+df
+
+
+
+# Create and plot the efficient frontier object
+ef <- efficientFrontierTwoAsset(R_C, R_S, sigma_C, sigma_S, rho)
+plot(ef)
+
+
+
+minRisk(R_C, R_S, sigma_C, sigma_S, rho)
+
+
+
+X_C <- seq(from=0, to=1, by=0.01)
+R_P <- portReturnTwoAsset(R_C, R_S, X_C)
+# rho = 1
+plot(portSDTwoAsset(R_C, R_S, X_C, sigma_C, sigma_S, rho=1), R_P,
+ main="Portfolio Return and Standard Deviation",
+ type="l", xlim=c(0, 0.08), ylim=c(0, 0.18),
+ xlab="Portfolio Standard Deviation",
+ ylab="Portfolio Expected Return", cex.lab=0.8)
+# rho = -1
+lines(portSDTwoAsset(R_C, R_S, X_C, sigma_C, sigma_S, rho=-1), R_P, col="blue")
+# rho = 0
+lines(portSDTwoAsset(R_C, R_S, X_C, sigma_C, sigma_S, rho=0), R_P, col="red")
+# rho = 0.5
+lines(portSDTwoAsset(R_C, R_S, X_C, sigma_C, sigma_S, rho=0.5), R_P, col="green")
+legend("topleft", legend=c(legend=expression(paste(rho, " = 1")),
+ legend=expression(paste(rho, " = -1")),
+ legend=expression(paste(rho, " = 0")),
+ legend=expression(paste(rho, " = 0.5"))),
+ col=c("black", "blue", "red", "green"),
+ lty=rep(1, 4), bty="n")
+
+
+
+# correlation coefficient
+rho <- 0.5
+
+# Fraction of the portfolio invested in Colonel Motors
+X_C <- seq(from=-1, to=2, by=0.2)
+
+# Calculate the portfolio expected return
+R_P <- portReturnTwoAsset(R_C, R_S, X_C)
+
+# Calculate the portfolio standard deviation
+sigma_P <- portSDTwoAsset(R_C, R_S, X_C, sigma_C, sigma_S, rho)
+
+# Combine the portfolio returns and standard deviations in a data.frame object.
+# Note that this is the transpose t() of the data.frame object to match the
+# layout of the tables in the FRM book.
+df <- t(data.frame(R_P=R_P, sigma_P=sigma_P))
+colnames(df) <- X_C
+print(df, digits=4)
+
+
+
+# Create and plot the efficient frontier object
+ef_short <- efficientFrontierTwoAsset(R_C, R_S, sigma_C, sigma_S, rho,
+ allowShorting=TRUE)
+ef_long <- efficientFrontierTwoAsset(R_C, R_S, sigma_C, sigma_S, rho)
+plot(ef_short)
+lines(ef_long$efficient_frontier[,2], ef_long$efficient_frontier[,1],
+ lty=2, lwd=2, col="blue")
+
+
+
+# Estimated inputs for equity
+R_SP <- 0.125
+sigma_SP <- 0.149
+
+# Estimated inputs for bonds
+R_B <- 0.06
+sigma_B <- 0.048
+
+# Estimated correlation between equity and bonds
+rho <- 0.45
+
+
+
+# Calculate the allocation and values for the minimum variance portfolio
+minRisk(R_SP, R_B, sigma_SP, sigma_B, rho)
+
+
+
+ef <- efficientFrontierTwoAsset(R_SP, R_B, sigma_SP, sigma_B, rho,
+ rf=0.05, allowShorting=TRUE)
+plot(ef)
+
+
+
+# Estimated inputs for domestic portfolio
+R_SP <- 0.125
+sigma_SP <- 0.149
+
+# Estimated inputs for international portfolio
+R_int <- 0.105
+sigma_int <- 0.14
+
+# Estimated correlation between domestic and international portfolio
+rho <- 0.33
+
+
+
+# Calculate the allocation and values for the minimum variance portfolio
+minRisk(R_SP, R_int, sigma_SP, sigma_int, rho)
+
+
+
+ef <- efficientFrontierTwoAsset(R_SP, R_int, sigma_SP, sigma_int, rho, rf=0.05)
+plot(ef)
+
+
+
+data(crsp_weekly)
+R_large <- largecap_weekly[,1:5]
+R_mid <- midcap_weekly[,1:5]
+R_small <- smallcap_weekly[,1:5]
+
+
+
+# Create and plot an efficient frontier of large cap stocks
+ef_large <- efficientFrontier(R_large)
+plot(ef_large, main="Large Cap Efficient Frontier", cexAssets=0.6)
+
+
+
+# Create an efficient frontier object with box constraints
+# Add box constraints such that each asset must have a weight greater than 15%
+# and less than 55%
+ef_large_box <- efficientFrontier(R_large, minBox=0.15, maxBox=0.55)
+plot(ef_large_box, main="Large Cap Efficient Frontier\nwith Box Constraints",
+ cexAssets=0.6)
+
+
+
+# Combine the large cap, mid cap, and small cap stocks
+R <- cbind(R_large, R_mid, R_small)
+
+# Specify the list of groups such that each market cap is in its own group
+groups <- list(c(1:5), c(6:10), c(11:15))
+
+# Specify the sum of weights for each group:
+# largecap: no more than 60% and no less than 30%
+# midcap: no more than 50% and no less than 20%
+# smallcap: no more than 40% and no less than 15%
+groupMin <- c(0.3, 0.2, 0.15)
+groupMax <- c(0.6, 0.5, 0.4)
+
+# Create the efficient frontier with the group constraints
+efGroup <- efficientFrontier(R, groupList=groups,
+ groupMin=groupMin,
+ groupMax=groupMax)
+plot(efGroup, labelAssets=FALSE)
+
+
Added: pkg/GARPFRM/demo/EstimatingVolatilitiesCorrelation.R
===================================================================
--- pkg/GARPFRM/demo/EstimatingVolatilitiesCorrelation.R (rev 0)
+++ pkg/GARPFRM/demo/EstimatingVolatilitiesCorrelation.R 2014-04-15 23:12:16 UTC (rev 167)
@@ -0,0 +1,138 @@
+
+
+library(knitr)
+opts_chunk$set(cache=TRUE, tidy=FALSE, warning=FALSE, fig.width=5, fig.height=5)
+
+
+
+suppressPackageStartupMessages(library(GARPFRM))
+data(crsp_weekly)
+
+# Use the weekly MSFT returns
+R <- largecap_weekly[, "MSFT"]
+
+
+
+lambda <- 0.94
+initialWindow <- 15
+volEst <- EWMA(R, lambda, initialWindow, type="volatility")
+volEst
+
+
+
+vol <- realizedVol(R, n=5)
+
+
+
+plot(vol, main="EWMA Volatility Estimate vs. Realized Volatility")
+lines(volEst$estimate, col="red")
+legend("topright", legend=c("Realized Volatility", "EWMA Volatility Estimate"),
+ col=c("black", "red"), lty=c(1,1), cex=0.8, bty="n")
+
+
+
+# Estimate lambda
+# Use initialWindow = 15 for the EWMA volatility estimate and
+# n = 5 to calculate the realized volatility
+lambda <- estimateLambdaVol(R, initialWindow, n=5)
+lambda
+
+volEst2 <- EWMA(R, lambda, initialWindow, type="volatility")
+volEst2
+
+
+
+# Realized volatility
+plot(vol, main="EWMA Volatility Estimate vs. Realized Volatility")
+# EWMA volatility estimate, lambda = 0.94
+lines(volEst$estimate, col="red")
+# EWMA volatility estimate, lambda = 0.0.7359253
+lines(volEst2$estimate, col="blue")
+legend("topright", legend=c("Realized Volatility",
+ "EWMA Volatility, lambda = 0.94",
+ "EWMA Volatility, lambda = 0.7359253"),
+ col=c("black", "red", "blue"), lty=c(1, 1, 1), cex=0.7, bty="n")
+
+
+
+# Use the first 2 columns of the large cap weekly returns
+R <- largecap_weekly[,1:2]
+initialWindow <- 52
+covEst <- EWMA(R, lambda=NULL, initialWindow, n=10, "covariance")
+covEst
+plot(covEst, main="EWMA Estimated Covariance")
+
+
+
+corEst <- EWMA(R, lambda=NULL, initialWindow, n=10, "correlation")
+corEst
+plot(corEst, main="EWMA Estimated Correlation")
+
+
+
+# Use the first 4 columns of the largecap_weekly dataset
+R <- largecap_weekly[,1:4]
+
+# calculate the sample covariance matrix
+sample_cov <- cov(R)
+sample_cov
+
+# EWMA covariance matrix estimate
+lambda <- 0.94
+initialWindow <- 52
+covEst <- EWMA(R, lambda, initialWindow, type="covariance")
+covEst
+
+
+
+# calculate the sample covariance matrix
+sample_cor <- cor(R)
+sample_cor
+
+# EWMA covariance matrix estimate
+lambda <- 0.94
+initialWindow <- 52
+corEst <- EWMA(R, lambda, initialWindow, type="correlation")
+corEst
+
+
+
+# Use the weekly MSFT returns
+R <- largecap_weekly[,"MSFT"]
+
+# Specify and fit the MSFT returns to a standard ARMA(0,0)-GARCH(1,1) model
+# Note that the default is ARMA(1,1)-GARCH(1,1) so we only need to change
+# the ARMA order. The default arguments were chosen to be consistent with the
+# default arguments in rugarch.
+model <- uvGARCH(R, armaOrder=c(0,0))
+
+# Get the fitted GARCH model
+fit <- getFit(model)
+
+# Get the coefficients
+coef(fit)
+
+# Show the summary results of the fit
+fit
+
+
+
+# n period ahead forecast
+forecast10 <- forecast(model, nAhead=10)
+forecast10
+
+
+
+plot(forecast10, which=3)
+
+
+
+model11 <- uvGARCH(R, armaOrder=c(0,0), outSample=100)
+forecast2 <- forecast(model11, nRoll=10)
+forecast2
+
+
+
+plot(forecast2, which=4)
+
+
Added: pkg/GARPFRM/demo/MonteCarloMethods.R
===================================================================
--- pkg/GARPFRM/demo/MonteCarloMethods.R (rev 0)
+++ pkg/GARPFRM/demo/MonteCarloMethods.R 2014-04-15 23:12:16 UTC (rev 167)
@@ -0,0 +1,187 @@
+
+
+library(knitr)
+opts_chunk$set(tidy=FALSE, warning=FALSE, fig.width=5, fig.height=5)
+
+
+
+suppressPackageStartupMessages(library(GARPFRM))
+# drift rate
+mu <- 0
+
+# volatility rate
+sigma <- 0.1
+
+# starting price
+S0 <- 100
+
+# number of steps
+N <- 100
+
+dt <- 1 / N
+
+# Generate N standard normal random variables
+set.seed(123)
+eps <- rnorm(N)
+
+# Allocate a vector to hold the prices
+S <- vector("numeric", N+1)
+S[1] <- S0
+
+# Precompute some of the terms
+mu_dt <- mu * dt
+sig_dt <- sigma * sqrt(dt)
+
+for(i in 2:length(S)){
+ S[i] <- S[i-1] + S[i-1] * (mu_dt + sig_dt * eps[i-1])
+}
+head(S)
+plot(S, main="Simulated Price Path", type="l")
+
+
+
+# Allocate a vector to hold the prices
+S1 <- vector("numeric", N+1)
+S1[1] <- S0
+
+# Precompute terms
+mu_sig_dt <- (mu - 0.5 * sigma^2) * dt
+sig_dt <- sigma * sqrt(dt)
+
+for(i in 2:length(S1)){
+ S1[i] <- S1[i-1] * exp(mu_sig_dt + sig_dt * eps[i-1])
+}
+head(S1)
+plot(S1, main="Simulated Price Path", type="l")
+
+
+
+mu <- 0.05
+sigma <- 0.15
+N <- 10000
+time <- 1
+steps <- 52
+startingValue <- 100
+
+# Run Monte Carlo simulation and store simulated price paths
+mcSim <- monteCarlo(mu, sigma, N, time, steps, startingValue)
+summary(endingPrices(mcSim))
+
+
+
+par(mfrow=c(2,1))
+plot(mcSim)
+plotEndingPrices(mcSim)
+par(mfrow=c(1,1))
+
+
+
+data(crsp_weekly)
+R.MSFT <- largecap_weekly[, "MSFT"]
+
+# Project number of periods ahead
+nAhead <- 5
+
+# Previous price
+S.p <- 25
+
+# Using a for loop
+bootS <- vector("numeric", nAhead)
+for(i in 1:nAhead){
+ bootS[i] <- S.p * (1 + sample(R.MSFT, 1, TRUE))
+ S.p <- bootS[i]
+}
+bootS
+
+# Vectorized solution
+S.p <- 25
+bootS1 <- S.p * cumprod(1 + sample(coredata(R.MSFT), nAhead, TRUE))
+bootS1
+
+
+
+# Number of boostrap replications
+rep <- 10000
+
+# Allocate vector to hold VaR statistic
+out <- vector("numeric", rep)
+for(i in 1:rep){
+ out[i] <- VaR(R.MSFT[sample.int(nrow(R.MSFT), replace=TRUE)],
+ method="historical")
+}
+
+# Bootstrapped VaR
+mean(out)
+
+# Standard error of Bootstrapped VaR
+sd(out)
+
+
+
+R <- largecap_weekly[,1:4]
+
+# function to calculate the annualized StdDev using the most recent n periods
+foo <- function(R, n){
+ StdDev.annualized(tail(R, n), geometric=TRUE)
+}
+
+bootFUN(R[,1], FUN="foo", n=104, replications=1000)
+
+# Bootstrap mean estimate.
+bootMean(R)
+
+# Bootstrap standard deviation estimate.
+bootSD(R)
+
+# Bootstrap standard deviation estimate using the StdDev function from
+# PerformanceAnalytics.
+bootStdDev(R)
+
+# Bootstrap simpleVolatility estimate.
+bootSimpleVolatility(R)
+
+# Bootstrap correlation estimate.
+bootCor(R)
+
+# Bootstrap covariance estimate.
+bootCov(R)
+
+# Bootstrap Value-at-Risk (VaR) estimate using the VaR function from
+# PerformanceAnalytics.
+bootVaR(R, p=0.9, method="historical", invert=FALSE)
+bootVaR(R, p=0.9, method="gaussian", invert=FALSE)
+
+# Bootstrap Expected Shortfall (ES) estimate using the ES function from
+# PerformanceAnalytics. Also known as Conditional Value-at-Risk (CVaR) and
+# Expected Tail Loss (ETL).
+bootES(R, p=0.9, method="historical")
+bootES(R, p=0.9, method="gaussian")
+
+
+
+# Register multicore parallel backend with 3 cores
+# Note that this example does not work on Windows
+# Windows users should use doSNOW
+# library(doMC)
+# registerDoMC(3)
+
+# Estimate VaR via bootstrap
+# bootVaR(R[,1], p=0.9, method="historical", replications=1000, parallel=TRUE)
+
+
+# Benchmark the performance of running the bootstrap in parallel
+# Bootstrap VaR with parallel=TRUE
+# bootPar <- function(){
+# bootVaR(R[,1], p=0.9, method="historical", replications=5000, parallel=TRUE)
+# }
+
+# Bootstrap VaR with parallel=FALSE
+# bootSeq <- function(){
+# bootVaR(R[,1], p=0.9, method="historical", replications=5000, parallel=FALSE)
+# }
+
+# Benchmark these functions
+# library(rbenchmark)
+# benchmark(bootPar(), bootSeq(), replications=1)[,1:4]
+
+
Added: pkg/GARPFRM/demo/PerformanceMeasures.R
===================================================================
--- pkg/GARPFRM/demo/PerformanceMeasures.R (rev 0)
+++ pkg/GARPFRM/demo/PerformanceMeasures.R 2014-04-15 23:12:16 UTC (rev 167)
@@ -0,0 +1,116 @@
+
+
+library(knitr)
+opts_chunk$set(cache=TRUE, tidy=FALSE, warning=FALSE, fig.width=5, fig.height=5)
+
+
+
+# Load the GARPFRM package and the CRSP dataset.
+suppressPackageStartupMessages(library(GARPFRM))
+data(crsp.short)
+
+# Market returns
+R.market <- largecap.ts[, "market"]
+
+# risk free rate
+rf <- largecap.ts[,"t90"]
+
+# The portfolio we will consider is an equal weight portfolio of the first
+# 10 assets in largecap.ts
+R.portfolio <- Return.portfolio(largecap.ts[,1:10])
+
+# Precompute excess returns
+R.Ex.portfolio <- R.portfolio - rf
+R.Ex.market <- R.market - rf
+
+
+
+# Compute portfolio beta using the covariance of the portfolio and benchmark
+# portfolio returns divided by the variance of the market portfolio returns
+cov(R.Ex.portfolio, R.Ex.market) / var(R.Ex.market)
+
+# Compute beta using CAPM
+fit <- CAPM(R.Ex.portfolio, R.Ex.market)
+getBetas(fit)
+
+# We can also directly use the CAPM.beta function from PerformanceAnalytics
+CAPM.beta(R.portfolio, R.market, rf)
+
+
+
+# Treynor ratio for portfolio and market
+
+# Treynor Ratio for portfolio
+TreynorRatio(R.portfolio, R.market, rf)
+
+# Treynor Ratio for market
+TreynorRatio(R.market, R.market, rf)
+
+
+
+# Compute Sharpe and annualized Sharpe Ratio
+# Sub-period Sharpe Ratio
+SharpeRatio(R.portfolio, rf, FUN="StdDev")
+
+# Annualized Sharpe Ratio
+SharpeRatio.annualized(R.portfolio, rf)
+
+
+
+SharpeRatio(R.portfolio, rf, p=0.95, FUN=c("VaR", "ES"))
+
+
+
+# Compute Jensen's alpha by carrying out a linear regression
+fit <- lm(R.Ex.portfolio ~ R.Ex.market)
+alpha <- coef(fit)[1]
+p_value <- coef(summary(fit))[1,4]
+summary(fit)
+
+# Compute Jensen's alpha with PerformanceAnalytics function
+CAPM.jensenAlpha(R.portfolio, R.market, mean(rf))
+
+# Replicate CAPM.jensenAlpha
+# Compute annualized returns
+R.P <- Return.annualized(R.portfolio)
+R.M <- Return.annualized(R.market)
+# Compute the CAPM beta
+beta <- CAPM.beta(R.portfolio, R.market, mean(rf))
+
+# Jensen's alpha
+R.P - mean(rf) - beta * (R.M - mean(rf))
+
+
+
+# Compute Tracking Error
+TrackingError(R.portfolio, R.market)
+
+# Replicate TrackingError
+sd(R.portfolio - R.market) * sqrt(12)
+
+
+
+# Compute Information Ratio
+# InformationRatio = ActivePremium / TrackingError
+# Active Premium = Investment's annualized return - Benchmark's annualized return
+InformationRatio(R.portfolio, R.market)
+
+# Replicate the Information Ratio computation
+activePremium <- Return.annualized(R.portfolio) - Return.annualized(R.market)
+trackingError <- TrackingError(R.portfolio, R.market)
+activePremium / trackingError
+
+
+
+
+# Compute Downside Deviation
+MAR <- 0
+# PA computation of Downside Deviation
+DownsideDeviation(R.portfolio, MAR)
+
+
+
+# Compute Sortino Ratio
+SortinoRatio(R.portfolio, MAR)
+
+
Added: pkg/GARPFRM/demo/QuantifyingVolatilityVaRModels.R
===================================================================
--- pkg/GARPFRM/demo/QuantifyingVolatilityVaRModels.R (rev 0)
+++ pkg/GARPFRM/demo/QuantifyingVolatilityVaRModels.R 2014-04-15 23:12:16 UTC (rev 167)
@@ -0,0 +1,300 @@
+
+
+library(knitr)
+opts_chunk$set(cache=TRUE, tidy=FALSE, warning=FALSE, fig.width=5, fig.height=5)
+
+
+
+# Load the package and data
+library(GARPFRM)
+data(crsp_weekly)
+R <- largecap_weekly
+R.COP <- R[,"COP"]
+
+
+
+plot(R.COP, main="COP Weekly Returns")
+
+
+
+hist(R.COP, breaks=50, main="Histogram of COP Returns",
+ col="lightblue", probability=TRUE)
+lines(density(R.COP), lwd=2)
+curve(dnorm(x, mean=mean(R.COP), sd=sd(R.COP)),
+ add=TRUE, col="red", lty=2, lwd=2)
+rug(R.COP)
+legend("topleft", legend=c("density", "normal"),
+ col=c("black", "red"), lty=c(1, 2), bty="n")
+
+
+
+hist(R.COP, breaks=50, main="Histogram of COP Returns",
+ col="lightblue", probability=TRUE,
+ xlim=c(min(density(R.COP)$x), -0.05))
+lines(density(R.COP), lwd=2)
+curve(dnorm(x, mean=mean(R.COP), sd=sd(R.COP)),
+ add=TRUE, col="red", lty=2, lwd=2)
+rug(R.COP)
+legend("topleft", legend=c("density", "normal"),
+ col=c("black", "red"), lty=c(1, 2), bty="n")
+
+
+
+chart.QQPlot(R.COP)
+
+
+
+# Compute rolling standard deviation estimate
+SD6 <- rollSD(R.COP, 6)
+SD13 <- rollSD(R.COP, 13)
+SD52 <- rollSD(R.COP, 52)
+# Plot rolling standard deviation estimates
+plot(SD6, type="n", main="Rolling Standard Deviation",
+ ylab="standard deviation")
+lines(SD6, col="blue", lty=2)
+lines(SD13, col="red")
+lines(SD52, lwd=2)
+legend("topleft", legend=c("rollSD (6)", "rollSD(13)", "rollSD(52)"),
+ bty="n", lty=c(2, 1, 1), col=c("blue", "red", "black"), cex=0.8)
+
+
+
+# Estimating volatility
+# EWMA Model
+initialWindow <- 100
+n <- 13
+type <- "volatility"
+
+# Fit an EWMA model to estimate volatility
+# Choose an optimal lambda parameter that minimizes the mean squared error
+# betwen realized volatility and the EWMA model volatility estimate
+ewmaModel <- EWMA(R.COP, lambda=NULL, initialWindow, n, type)
+ewmaModel
+
+# One period ahead forecast
+ewmaVolForecast <- forecast(ewmaModel)
+
+# VaR using EWMA volatility forecast
+# Here we assume the expected value of returns is simply the mean of returns
+mean(R.COP) + as.numeric(ewmaVolForecast) * qnorm(0.05)
+
+
+
+# Specify and fit a GARCH Model
+garchModel <- uvGARCH(R.COP, armaOrder=c(0,0))
+
+# One period ahead forecast of GARCH model
+garchForecast <- forecast(garchModel, 1)
+
+# VaR forecast using GARCH Model
+fitted(garchForecast) + sigma(garchForecast) * qnorm(0.05)
+
+
+
+# Historical VaR estimate at the 5% level
+historicalVaR5 <- VaR(R.COP, p=0.95, method="historical")
+
+# VaR estimate assuming a normal distribution
+normalVaR5 <- VaR(R.COP, p=0.95, method="gaussian")
+
+# Bootstrapped historical VaR estimate
+bootHistVaR5 <- bootVaR(R.COP, p=0.95, method="historical")
+
+rnames <- c("Historical", "Normal", "Bootstrap Historical")
+matrix(c(historicalVaR5, normalVaR5, bootHistVaR5[1,]),
+ nrow=3, ncol=1, dimnames=list(rnames, "VaR (5%)"))
+
+
+
+hist(R.COP, main="Histogram of COP returns", breaks=50,
+ col="blue", probability=TRUE)
+lines(density(R.COP), lwd=2)
+curve(dnorm(x, mean=mean(R.COP), sd=sd(R.COP)),
+ add=TRUE, col="red", lty=2, lwd=2)
+rug(R.COP)
+arrows(historicalVaR5, 0, historicalVaR5, 6, code=1, lwd=2)
+text(historicalVaR5, 6, labels="historical VaR (5%)", pos=2, cex=0.7)
+
+
+
+# Estimating VaR at the 1% level
+historicalVaR1 <- VaR(R.COP, p=0.99, method="historical")
+normalVaR1 <- VaR(R.COP, p=0.99, method="gaussian")
+bootHistVaR1 <- bootVaR(R.COP, p=0.99, method="historical")
+
+matrix(c(historicalVaR1, normalVaR1, bootHistVaR1[1,]),
+ nrow=3, ncol=1, dimnames=list(rnames, "VaR (1%)"))
+
+
+
+hist(R.COP, main="Histogram of COP returns", breaks=50,
+ col="blue", probability=TRUE)
+lines(density(R.COP), lwd=2)
+curve(dnorm(x, mean=mean(R.COP), sd=sd(R.COP)),
+ add=TRUE, col="red", lty=2, lwd=2)
+rug(R.COP)
+arrows(historicalVaR1, 0, historicalVaR1, 4, code=1, lwd=2)
+text(historicalVaR1, 4, labels="Historical VaR (1%)", pos=2, cex=0.7)
+arrows(normalVaR1, 0, normalVaR1, 6, code=1, lwd=2)
+text(normalVaR1, 6, labels="Normal VaR (1%)", pos=2, cex=0.7)
+
+
+
+# Asset returns
+R <- largecap_weekly[, 1:10]
+
+# Lookback period
+K <- 52
+
+# Set the weights K periods ago
+weights <- xts(matrix(rep(1 / 10, 10), nrow=1), index(R)[nrow(R) - K])
+
+# Calculate the portfolio returns for the most recent K periods
+R.portfolio <- Return.rebalancing(R, weights)
+
+
+
+par(mfrow=c(2,1))
+# Portfolio returns
+plot(R.portfolio, main="Portfolio Returns")
+
+# Histogram of portfolio returns
+hist(R.portfolio, main="Histogram of Portfolio returns", breaks=50,
+ col="blue", probability=TRUE)
+lines(density(R.portfolio), lwd=2)
+curve(dnorm(x, mean=mean(R.portfolio), sd=sd(R.portfolio)),
+ add=TRUE, col="red", lty=2, lwd=2)
+rug(R.portfolio)
+par(mfrow=c(1,1))
+
+
+
+# Estimate the VaR of the portfolio using the last 52 periods
+# Historical VaR estimate
+portfVaR.HS <- VaR(R.portfolio, p=0.95, method="historical")
+
+# Bootstrapped VaR estimate
+portfVaR.BootHS <- bootVaR(R.portfolio, p=0.95, method="historical")
+
+# Normal VaR estimate
+portfVaR.normal <- VaR(R.portfolio, p=0.95, method="gaussian")
+
+
+
+# Use the most recent 52 periods to compute the sample covariance matrix
+sampleCov <- cov(tail(R, 52))
+# Convert the xts object of weights to a matrix
+weights <- matrix(weights, ncol=1)
+
+# Compute the portfolio VaR estimate using the VarCov approach with sample
+# covariance matrix
+portfVaR.cov <- sqrt(t(weights) %*% sampleCov %*% weights) * qnorm(0.05)
+
+# Use EWMA model to compute variance covariance matrix
+# EWMA model to compute variance covariance matrix
+ewmaCov <- EWMA(tail(R, 52), lambda=0.9, initialWindow=10,
+ type="covariance")$estimate
+
+# Compute the portfolio VaR estimate using the VarCov approach with EWMA
+# model estimated covariance matrix
+portfVaR.ewmaCov <- sqrt(t(weights) %*% ewmaCov %*% weights) * qnorm(0.05)
+
+
+
+# Compute rolling correlation estimates
+cor13 <- rollCor(R[,1:2], 13)
+cor26 <- rollCor(R[,1:2], 26)
+
+# Compute rolling covariance estimates
+cov13 <- rollCov(R[,1:2], 13)
+cov26 <- rollCov(R[,1:2], 26)
+
+# Plot rolling correlation estimates
+plot(cor13, type="n", main="Rolling Correlation",
+ ylab="correlation")
+lines(cor13, col="blue")
+lines(cor26, col="red")
+legend("topleft", legend=c("rollCor(13)", "rollCor(26)"),
+ bty="n", lty=c(1, 1), col=c("blue", "red"), cex=0.8)
+
+# Plot rolling covariance estimates
+plot(cov13, type="n", main="Rolling Covariance",
+ ylab="covariance")
+lines(cov13, col="blue")
+lines(cov26, col="red")
+legend("topleft", legend=c("rollCov(13)", "rollCov(26)"),
+ bty="n", lty=c(1, 1), col=c("blue", "red"), cex=0.8)
+
+
+
+# Portfolio VaR estimate with EWMA model
+ewmaModel <- EWMA(R.portfolio, lambda=NULL, initialWindow=10, n, type)
+# One period ahead forecast
+ewmaVolForecast <- forecast(ewmaModel)
+
+# VaR using EWMA volatility forecast
+# Here we assume the expected value of returns is simply the mean of returns
+portfVaR.EWMA <- mean(R.portfolio) + as.numeric(ewmaVolForecast) * qnorm(0.05)
+
+
+
+dfVaR <- t(data.frame(portfVaR.HS, portfVaR.BootHS[1,], portfVaR.normal,
+ portfVaR.cov, portfVaR.ewmaCov, portfVaR.EWMA))
+rownames(dfVaR) <- c("portfVaR.HS", "portfVaR.BootHS", "portfVaR.normal",
+ "portfVaR.cov", "portfVaR.ewmaCov", "portfVaR.EWMA")
+dfVaR
+
+
+
+R <- largecap_weekly[, 1:10]
+
+# Annual rebalance dates
+rebalanceDates <- index(R)[endpoints(index(R), on="years")]
+
+# Create an xts object of weights at the specified rebalance dates
+weights <- xts(matrix(1 / 10, nrow=length(rebalanceDates),
+ ncol=10), rebalanceDates)
+
+# Calculate the aggregate portfolio return
+R.portfolio <- Return.rebalancing(R, weights)
+
+
+
+par(mfrow=c(2,1))
+# Portfolio returns
+plot(R.portfolio, main="Portfolio Returns")
+
+# Histogram of portfolio returns
+hist(R.portfolio, main="Histogram of Portfolio returns", breaks=50,
+ col="blue", probability=TRUE)
+lines(density(R.portfolio), lwd=2)
+curve(dnorm(x, mean=mean(R.portfolio), sd=sd(R.portfolio)),
+ add=TRUE, col="red", lty=2, lwd=2)
+rug(R.portfolio)
+par(mfrow=c(1,1))
+
+
+
+garchModel <- uvGARCH(R.portfolio, armaOrder=c(0,0))
+btVaR.GARCH <- backtestVaR.GARCH(garchModel, p=0.95, refitEvery=5, window=100)
+btVaR.GARCH
+
+
+
+# Plot the GARCH VaR backtest
+plot(btVaR.GARCH, pch=20, legendLoc="topright")
+
+
+
+# Run a VaR backtest on portfolio returns
+# Compute VaR estimate using gaussian, historical, and modified methods
+backtest <- backtestVaR(R.portfolio, window=100, p=0.95,
+ method=c("gaussian", "historical", "modified"))
+backtest
+
+
+
+# plot the VaR backtest
+plot(backtest, pch=18, legendLoc="topright")
+
+
Added: pkg/GARPFRM/demo/QuantitativeAnalysisBasics.R
===================================================================
--- pkg/GARPFRM/demo/QuantitativeAnalysisBasics.R (rev 0)
+++ pkg/GARPFRM/demo/QuantitativeAnalysisBasics.R 2014-04-15 23:12:16 UTC (rev 167)
@@ -0,0 +1,350 @@
+
+
+suppressMessages(library(GARPFRM))
+suppressMessages(library(lattice))
+suppressMessages(library(pcaPP))
+suppressMessages(library(hexbin))
+data(crsp.short)
+
+# Use the first 6 stocks in largecap.ts
+crsp.returns <- largecap.ts[, 1:6]
+
+# Get the names of the stocks and view the first 5 rows of crsp.returns
+colnames(crsp.returns)
+head(crsp.returns, 5)
+
+
+
+xyplot(crsp.returns, scale=list(y="same"), main="Monthly Returns")
+
+
+
+boxplot(coredata(crsp.returns), pch=20, main="Monthly Returns")
+
+
+
+# Extract the column labeled "market" to get the market returns
+MKT.ret <- largecap.ts[, "market"]
+
+
+
+plot(MKT.ret, main="Market Monthly Returns")
+
+
+
+hist(MKT.ret, probability=TRUE, main="Histogram of Market Returns",
+ col="lightblue", ylim=c(0, 10))
+lines(density(MKT.ret), lty=2)
+rug(MKT.ret)
+legend("topleft", legend="kernel density estimate", lty=2,
+ cex=0.8, bty="n")
+
+
+
+# Plot the kernel density estimate of market returns
+plot(density(MKT.ret), main="Density of Market Returns")
+rug(MKT.ret)
+# sample estimates
+curve(dnorm(x, mean=mean(MKT.ret), sd=sd(MKT.ret)),
+ add=TRUE, col="red", lty=2, lwd=2)
+# robust estimates
+curve(dnorm(x, mean=median(MKT.ret), sd=mad(MKT.ret)),
+ add=TRUE, col="blue", lty=2, lwd=2)
+legend("topleft", legend=c("kernel density estimate", "normal density estimate",
+ "robust normal density estimate"),
+ col=c("black", "red", "blue"), lty=c(1, 2, 2), bty="n", cex=0.8)
+
+
+
+chart.QQPlot(MKT.ret, envelope=0.95, pch=18, main="Market Returns QQ Plot",
+ xlab="Theoretical Normal Quantiles")
+legend("topleft", legend=c("Quartile-Pairs Line", "95% Confidence Envelope"),
+ col=c("blue", "blue"), lty=c(1, 2), cex=0.8, bty="n")
+
+
+
+shapiro.test(coredata(MKT.ret))
+
+
+
+# Sample mean of SPY return
+mean(MKT.ret)
+
+# Sample Variance of SPY returns
+var(MKT.ret)
+
+# Sample standard deviation of SPY returns
+sd(MKT.ret)
+
+# Standard error of SPY returns
+sd(MKT.ret) / sqrt(nrow(MKT.ret))
+
+# Sample skewness of SPY returns.
+# See ?skewness for additional methods for calculating skewness
+skewness(MKT.ret, method="sample")
+
+# Sample kurtosis of SPY returns.
+# See ?kurtosis for additional methods for calculating kurtosis
+kurtosis(MKT.ret, method="sample")
+
+# Summary statistics of SPY returns
+summary(MKT.ret)
+
+# Sample quantiles of SPY returns
+quantile(MKT.ret, probs=c(0, 0.25, 0.5, 0.75, 1))
+
+
+
+pairs(coredata(crsp.returns), pch=20, col=rgb(0,0,100,50,maxColorValue=255))
+hexplom(coredata(crsp.returns), varname.cex=0.75)
+
+
+
+# Sample correlation of returns
+cor(crsp.returns)
+
+
+
+suppressWarnings(plotcov(cor(crsp.returns),
+ method1="Sample Correlation Estimate"))
+
+
+
+# Sample covariance of returns
+cov(crsp.returns)
+
+
+
+suppressWarnings(plotcov(cov(crsp.returns),
+ method1="Sample Covariance Estimate"))
+
+
+
+curve(dnorm(x), from=-4, to=4, main="Standard Normal pdf")
+
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/uwgarp -r 167
More information about the Uwgarp-commits
mailing list