[Returnanalytics-commits] r3677 - in pkg/Meucci: . R demo man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Jun 16 17:39:42 CEST 2015
Author: xavierv
Date: 2015-06-16 17:39:42 +0200 (Tue, 16 Jun 2015)
New Revision: 3677
Modified:
pkg/Meucci/DESCRIPTION
pkg/Meucci/R/ButterflyTradingFunctions.R
pkg/Meucci/R/EntropyProg.R
pkg/Meucci/demo/ButterflyTrading.R
pkg/Meucci/demo/Prior2Posterior.R
pkg/Meucci/man/ComputeCVaR.Rd
pkg/Meucci/man/HorizonPricing.Rd
pkg/Meucci/man/LongShortMeanCVaRFrontier.Rd
pkg/Meucci/man/ViewCurveSlope.Rd
pkg/Meucci/man/ViewImpliedVol.Rd
pkg/Meucci/man/ViewRealizedVol.Rd
Log:
Fixed and formatted AnalyticalvsNumerical and ButterflyTrading codes and its dependencies (including EntropyProg)
Modified: pkg/Meucci/DESCRIPTION
===================================================================
--- pkg/Meucci/DESCRIPTION 2015-06-16 14:22:14 UTC (rev 3676)
+++ pkg/Meucci/DESCRIPTION 2015-06-16 15:39:42 UTC (rev 3677)
@@ -37,8 +37,8 @@
quadprog,
kernlab,
nloptr,
+ limSolve,
Suggests:
- limSolve,
Matrix,
MASS,
reshape2,
Modified: pkg/Meucci/R/ButterflyTradingFunctions.R
===================================================================
--- pkg/Meucci/R/ButterflyTradingFunctions.R 2015-06-16 14:22:14 UTC (rev 3676)
+++ pkg/Meucci/R/ButterflyTradingFunctions.R 2015-06-16 15:39:42 UTC (rev 3677)
@@ -1,19 +1,22 @@
# In order of appearance in the demo script ButterflyTrading.R
-MapVol = function( sig , y , K , T )
-{
- # in real life a and b below should be calibrated to security-specific time series
-
- a = -0.00000000001
- b = 0.00000000001
-
- s = sig + a/sqrt(T) * ( log(K) - log(y) ) + b/T*( log(K) - log(y) )^2
-
- return( s )
+MapVol <- function(sig, y, K, T) {
+ # in real life a and b below should be calibrated to security-specific time
+ # series
+
+ a <- -0.00000000001
+ b <- 0.00000000001
+
+ s <- sig + a / sqrt(T) * (log(K) - log(y)) + b / T * (log(K) - log(y)) ^ 2
+
+ return(s)
}
-#' Compute the pricing in the horizon as it appears in A. Meucci, "Fully Flexible Views: Theory and Practice",
-#' The Risk Magazine, October 2008, p 100-106.
+#' @title Compute the pricing in the horizon.
+#'
+#' @description Compute the pricing in the horizon as it appears in A. Meucci,
+#' "Fully Flexible Views: Theory and Practice", The Risk Magazine, October
+#' 2008, p 100-106.
#'
#' @param Butterflies List of securities with some analytics computed.
#' @param X Panel of joint factors realizations
@@ -21,125 +24,153 @@
#' @return PnL Matrix of profit and loss scenarios
#'
#' @references
-#' A. Meucci, "Fully Flexible Views: Theory and Practice" \url{http://www.symmys.com/node/158}
-#' See Meucci script for "ButterflyTrading/HorizonPricing.m"
+#' A. Meucci, "Fully Flexible Views: Theory and Practice"
+#' \url{http://www.symmys.com/node/158}, see Meucci script for
+#' "ButterflyTrading/HorizonPricing.m"
#'
-#' @author Ram Ahluwalia \email{ram@@wingedfootcapital.com} and Xavier Valls \email{flamejat@@gmail.com}
+#' @author Ram Ahluwalia \email{ram@@wingedfootcapital.com} and Xavier Valls
+#' \email{xaviervallspla@@gmail.com}
#' @export
-HorizonPricing = function( Butterflies , X )
-{
- r = 0.04 # risk-free rate
- tau = 1/252 # investment horizon
-
- # factors: 1. 'MSFT_close' 2. 'MSFT_vol_30' 3. 'MSFT_vol_91' 4. 'MSFT_vol_182'
- # securities: 1. 'MSFT_vol_30' 2. 'MSFT_vol_91' 3. 'MSFT_vol_182'
-
+HorizonPricing <- function(Butterflies, X) {
+ r <- 0.04 # risk-free rate
+ tau <- 1 / 252 # investment horizon
+
+ # factors: 1. 'MSFT_close' 2. 'MSFT_vol_30'
+ # 3. 'MSFT_vol_91' 4. 'MSFT_vol_182'
+ # securities: 1. 'MSFT_vol_30' 2. 'MSFT_vol_91'
+ # 3. 'MSFT_vol_182'
+
# create a new row called DlnY and Dsig
- # create a new row called 'DlnY'. Assign the first row (vector) of X to this DlnY for the 1:3 securities
- for ( s in 1:3 ) { Butterflies[[s]]$DlnY = X[ , 1 ] }
-
+ # create a new row called 'DlnY'. Assign the first row (vector) of X to this
+ # DlnY for the 1:3 securities
+ for (s in 1:3) {
+ Butterflies[[s]]$DlnY <- X[, 1]
+ }
+
# assign the 2nd row of X to a new element called Dsig
- Butterflies[[1]]$Dsig=X[ , 2 ]
- Butterflies[[2]]$Dsig=X[ , 3 ]
- Butterflies[[3]]$Dsig=X[ , 4 ]
-
- # factors: 5. 'YHOO_close' 6. 'YHOO_vol_30' 7. 'YHOO_vol_91' 8. 'YHOO_vol_182'
- # securities: 4. 'YHOO_vol_30' 5. 'YHOO_vol_91' 6. 'YHOO_vol_182'
- for ( s in 4:6 ) { Butterflies[[s]]$DlnY=X[ , 5 ] }
-
- Butterflies[[4]]$Dsig=X[ , 6 ]
- Butterflies[[5]]$Dsig=X[ , 7 ]
- Butterflies[[6]]$Dsig=X[ , 8 ]
-
- # factors: # 9. 'GOOG_close' 10. 'GOOG_vol_30' 11. 'GOOG_vol_91' 12. 'GOOG_vol_182'
- # securities: 7. 'GOOG_vol_30' 8. 'GOOG_vol_91' 9. 'GOOG_vol_182'
- for ( s in 7:9 ) { Butterflies[[s]]$DlnY=X[ , 9 ] }
-
- Butterflies[[7]]$Dsig=X[ , 10 ]
- Butterflies[[8]]$Dsig=X[ , 11 ]
- Butterflies[[9]]$Dsig=X[ , 12 ]
-
- PnL = matrix( NA , nrow = nrow(X) )
-
- for ( s in 1:length(Butterflies) )
- {
- Y = Butterflies[[s]]$Y_0 * exp(Butterflies[[s]]$DlnY)
- ATMsig = apply( cbind( Butterflies[[s]]$sig_0 + Butterflies[[s]]$Dsig , 10^-6 ) , 1 , max )
- t = Butterflies[[s]]$T - tau
- K = Butterflies[[s]]$K
- sig = MapVol(ATMsig , Y , K , t )
-
- ############# Ram's Code: Substituted with package's own functions #################################
+ Butterflies[[1]]$Dsig <- X[, 2]
+ Butterflies[[2]]$Dsig <- X[, 3]
+ Butterflies[[3]]$Dsig <- X[, 4]
+
+ # factors: 5. 'YHOO_close' 6. 'YHOO_vol_30' 7. 'YHOO_vol_91'
+ # 8. 'YHOO_vol_182'
+ # securities: 4. 'YHOO_vol_30' 5. 'YHOO_vol_91' 6. 'YHOO_vol_182'
+ for (s in 4:6) {
+ Butterflies[[s]]$DlnY <- X[, 5]
+ }
+
+ Butterflies[[4]]$Dsig <- X[, 6]
+ Butterflies[[5]]$Dsig <- X[, 7]
+ Butterflies[[6]]$Dsig <- X[, 8]
+
+ # factors: 9. 'GOOG_close' 10. 'GOOG_vol_30' 11. 'GOOG_vol_91'
+ # 12. 'GOOG_vol_182'
+ # securities: 7. 'GOOG_vol_30' 8. 'GOOG_vol_91' 9. 'GOOG_vol_182'
+ for (s in 7:9) {
+ Butterflies[[s]]$DlnY <- X[, 9]
+ }
+
+ Butterflies[[7]]$Dsig <- X[, 10]
+ Butterflies[[8]]$Dsig <- X[, 11]
+ Butterflies[[9]]$Dsig <- X[, 12]
+
+ PnL <- matrix(NA, nrow = nrow(X))
+
+ for (s in 1:length(Butterflies)) {
+ Y <- Butterflies[[s]]$Y_0 * exp(Butterflies[[s]]$DlnY)
+ ATMsig <- apply(cbind(Butterflies[[s]]$sig_0 + Butterflies[[s]]$Dsig,
+ 10 ^ -6), 1, max)
+ t <- Butterflies[[s]]$T - tau
+ K <- Butterflies[[s]]$K
+ sig <- MapVol(ATMsig, Y, K, t)
+
+ ############# Ram's Code: Substituted with package's own functions #########
#
- ## library(RQuantLib) # this function can only operate on one option at a time, so we use fOptions
- ##C = EuropeanOption( type = "call" , underlying = Y , strike = K , dividendYield = 0 , riskFreeRate = r , maturity = t , volatility = sig )$value
- ## P = EuropeanOption( type = "put" , underlying = Y , strike = K , dividendYield = 0 , riskFreeRate = r , maturity = t , volatility = sig )$value
-
+ ## library(RQuantLib) # this function can only operate on one option at a
+ ## time, so we use fOptions
+ ## C = EuropeanOption(type = "call", underlying = Y, strike = K,
+ ## dividendYield = 0, riskFreeRate = r, maturity = t,
+ ## volatility = sig)$value
+ ## P = EuropeanOption(type = "put", underlying = Y, strike = K,
+ ## dividendYield = 0, riskFreeRate = r, maturity = t,
+ ## volatility = sig)$value
## use fOptions to value options
- #library( fOptions )
- #C = GBSOption( TypeFlag = "c" , S = Y , X = K , r = r , b = 0 , Time = t , sigma = sig )
- #P = GBSOption( TypeFlag = "p" , S = Y , X = K , r = r , b = 0 , Time = t , sigma = sig )
- #
- ####################################################################################################
+ #library(fOptions)
+ #C = GBSOption(TypeFlag = "c", S = Y, X = K, r = r, b = 0, Time = t,
+ # sigma = sig)
+ #P = GBSOption(TypeFlag = "p", S = Y, X = K, r = r, b = 0, Time = t,
+ # sigma = sig)
+ ############################################################################
- BS = BlackScholesCallPutPrice( Y, K, r, sig, t )
-
- Butterflies[[s]]$P_T = BS$call + BS$put
- PnL = cbind( PnL , Butterflies[[s]]$P_T )
+ BS <- BlackScholesCallPutPrice(Y, K, r, sig, t)
+ Butterflies[[s]]$P_T <- BS$call + BS$put
+ PnL <- cbind(PnL, Butterflies[[s]]$P_T)
}
- PnL = PnL[ , -1 ]
-
- return( PnL )
+ PnL <- PnL[, -1]
+
+ return(PnL)
}
-ViewCurveSlopeTest = function( X , p )
-{
- J = nrow( X ) ; K = ncol( X )
-
+ViewCurveSlopeTest <- function(X, p) {
+ J <- nrow(X)
+
# constrain probabilities to sum to one...
- Aeq = matrix( 1, 1 , J )
- beq = matrix( 1 , nrow = 1 , ncol = 1 )
+ Aeq <- matrix(1, 1, J)
+ beq <- matrix(1, nrow = 1, ncol = 1)
browser()
# ...constrain the expectation...
- V = matrix( , nrow = nrow( X ) , ncol = 0 )
+ V <- matrix(, nrow = nrow(X), ncol = 0)
# Add 3 equality views
- V = cbind( V , X[ , 14 ] - X[ , 13 ] ) # View 1: spread on treasuries
- V = cbind( V , X[ , 14 ] - X[ , 13 ] ) # View 2: identical view (spread on treasuries)
- V = cbind( V , X[ , 6 ] - X[ , 5 ] ) # View 3: difference in YHOO Vol
- v = matrix( c( .0005 , 0 ) , nrow = ncol( V ) , ncol = 1 )
-
- Aeq = rbind( Aeq , t(V) )
-
- beq = rbind( beq , v )
-
+ V <- cbind(V, X[, 14] - X[, 13]) # View 1: spread on treasuries
+ V <- cbind(V, X[, 14] - X[, 13]) # View 2: identical view(spread on treasuries
+ V <- cbind(V, X[, 6] - X[, 5]) # View 3: difference in YHOO Vol
+ v <- matrix(c(.0005, 0), nrow = ncol(V), ncol = 1)
+
+ Aeq <- rbind(Aeq, t(V))
+
+ beq <- rbind(beq, v)
+
# add an inequality view
- # ...constrain the median...
- V = abs( X[ , 1 ] ) # absolute value of the log of changes in MSFT close prices (definition of realized volatility)
- V_Sort = sort( V , decreasing = FALSE ) # sorting of the abs value of log changes in prices from smallest to largest
- I_Sort = order( V )
-
- F = cumsum( p[ I_Sort ] ) # represents the cumulative sum of probabilities from ~0 to 1
-
- I_Reference = max( matlab:::find( F <= 3/5 ) ) # finds the (max) index corresponding to element with value <= 3/5 along the empirical cumulative density function for the abs log-changes in price
- V_Reference = V_Sort[ I_Reference ] # returns the corresponding abs log of change in price at the 3/5 of the cumulative density function
+ # ...constrain the median...
+ # absolute value of the log of changes in MSFT close prices (definition of
+ # realized volatility)
+ V <- abs(X[, 1])
+ # sorting of the abs value of log changes in prices from smallest to largest
+ V_Sort <- sort(V, decreasing = FALSE)
+ I_Sort <- order(V)
- I_Select = find( V <= V_Reference ) # finds all indices with value of abs log-change in price less than the reference value
- a = zeros( 1 , J )
- a[ I_Select ] = 1 # select those cases where the abs log-change in price is less than the 3/5 of the empirical cumulative density...
-
- A = a
- b = 0.5 # ... and assign the probability of these cases occuring as 50%. This moves the media of the distribution
-
+ # represents the cumulative sum of probabilities from ~0 to 1
+ F <- cumsum(p[I_Sort])
+
+ # finds the (max) index corresponding to element with value <= 3/5 along the
+ # empirical cumulative density function for the abs log-changes in price
+ I_Reference <- max(matlab:::find(F <= 3 / 5))
+ # returns the corresponding abs log of change in price at the 3/5 of the cdf
+ V_Reference <- V_Sort[I_Reference]
+
+ # finds all indices with value of abs log-change in price < reference value
+ I_Select <- find(V <= V_Reference)
+ a <- zeros(1, J)
+ # select those cases where the abs log-change in price is less than the 3/5 of
+ # the empirical cumulative density...
+ a[I_Select] <- 1
+
+ A <- a
+ # ... and assign the probability of these cases occuring as 50%.
+ # This moves the media of the distribution
+ b <- 0.5
+
# ...compute posterior probabilities
- p_ = EntropyProg( p , A , b , Aeq ,beq )
- return( p_ )
+ p_ <- EntropyProg(p, A, b, Aeq,beq)$p_
+ return(p_)
}
-
-#' Process the inequality view, as it appears in A. Meucci, "Fully Flexible Views: Theory and Practice",
-#' The Risk Magazine, October 2008, p 100-106.
+#' @title Process the inequality view
+#' @description Process the inequality view, as it appears in A. Meucci, "Fully
+#' Flexible Views: Theory and Practice", The Risk Magazine, October 2008,
+#' p 100-106.
#'
#' @param X : panel of joint factors realizations
#' @param p : vector of probabilities
@@ -147,40 +178,45 @@
#' @return p_ : vector of posterior probabilities
#'
#' @references
-#' A. Meucci, "Fully Flexible Views: Theory and Practice" \url{http://www.symmys.com/node/158}
+#' A. Meucci, "Fully Flexible Views: Theory and Practice"
+#' \url{http://www.symmys.com/node/158}
+#'
#' See Meucci script for "ButterflyTrading/ViewRealizedVol.m"
#'
-#' @author Ram Ahluwalia \email{ram@@wingedfootcapital.com} and Xavier Valls \email{flamejat@@gmail.com}
+#' @author Ram Ahluwalia \email{ram@@wingedfootcapital.com} and
+#' Xavier Valls \email{xaviervallspla@@gmail.com}
+#'
#' @export
-ViewImpliedVol = function( X , p )
-{
- # View 1 (inequality view): bearish on on 2m-6m implied volaility spread for Google
-
- J = nrow( X ) ;
- K = ncol( X );
-
+ViewImpliedVol <- function(X, p) {
+ # View 1 (inequality view): bearish on on 2m-6m implied volaility spread
+ #for Google
+
+ J <- nrow(X)
+
# constrain probabilities to sum to one...
- Aeq = matrix( 1, 1 , J )
- beq = 1
-
+ Aeq <- matrix(1, 1, J)
+ beq <- c(1)
+
# ...constrain the expectation...
- V = X[ , 12 ] - X[ , 11 ] # GOOG_vol_182 (6m implied vol) - GOOG_vol_91 (2m implied vol)
- m = mean( V )
- s = std( V )
-
- A = t( V )
- b = m - s
-
+ # GOOG_vol_182 (6m implied vol) - GOOG_vol_91 (2m implied vol)
+ V <- X[, 12] - X[, 11]
+ m <- mean(V)
+ s <- std(V)
+
+ A <- t(V)
+ b <- c(m - s)
+
# ...compute posterior probabilities
- p_ = EntropyProg( p , A , b , Aeq , beq )$p_
-
- return( p_ )
+ p_ <- EntropyProg(p, A, b, Aeq, beq)$p_
+
+ return(p_)
}
-#' Process the relative inequality view on median, as it appears in A. Meucci,
-#' "Fully Flexible Views: Theory and Practice", The Risk Magazine, October 2008,
-#' p 100-106
+#' @title Process the relative inequality view on median
+#' @description Process the relative inequality view on median, as it appears in
+#' A. Meucci, "Fully Flexible Views: Theory and Practice", The Risk Magazine,
+#' October 2008, p 100-106
#'
#' @param X : panel of joint factors realizations
#' @param p : vector of probabilities
@@ -188,53 +224,69 @@
#' @return p_ : vector of posterior probabilities
#'
#' @references
-#' A. Meucci, "Fully Flexible Views: Theory and Practice" \url{http://www.symmys.com/node/158}
+#' A. Meucci, "Fully Flexible Views: Theory and Practice"
+#' \url{http://www.symmys.com/node/158}
#' See Meucci script for "ButterflyTrading/ViewRealizedVol.m"
#'
-#' @author Ram Ahluwalia \email{ram@@wingedfootcapital.com} and Xavier Valls \email{flamejat@@gmail.com}
+#' @author Ram Ahluwalia \email{ram@@wingedfootcapital.com} and Xavier Valls
+#' \email{xaviervallspla@@gmail.com}
#' @export
-ViewRealizedVol = function( X , p )
-{
- # view 2 bullish on realized volatility of MSFT (i.e. absolute log-change in the underlying).
- # This is the variable such that, if larger than a threshold, a long position in the butterfly turns into a profit (e.g. Rachev 2003)
- # we issue a relative statement on the media comparing it with the third quintile implied by the reference market model
-
- library( matlab )
- J = nrow( X ) ; K = ncol( X )
-
+ViewRealizedVol <- function(X, p) {
+ # view 2 bullish on realized volatility of MSFT (i.e. absolute log-change in
+ # the underlying).
+ # This is the variable such that, if larger than a threshold, a long position
+ # in the butterfly turns into a profit (e.g. Rachev 2003)
+ # we issue a relative statement on the media comparing it with the third
+ # quintile implied by the reference market model
+
+ J <- nrow(X)
+
# constrain probabilities to sum to one...
- Aeq = matrix( 1, 1 , J )
- beq = 1
-
- # ...constrain the median...
- V = abs( X[ , 1 ] ) # absolute value of the log of changes in MSFT close prices (definition of realized volatility)
-
- V_Sort = sort( V , decreasing = FALSE ) # sorting of the abs value of log changes in prices from smallest to largest
- I_Sort = order( V )
-
- F = cumsum( p[ I_Sort ] ) # represents the cumulative sum of probabilities from ~0 to 1
-
- I_Reference = max( matlab:::find( F <= 3/5 ) ) # finds the (max) index corresponding to element with value <= 3/5 along the empirical cumulative density function for the abs log-changes in price
- V_Reference = V_Sort[ I_Reference ] # returns the corresponding abs log of change in price at the 3/5 of the cumulative density function
-
- I_Select = find( V <= V_Reference ) # finds all indices with value of abs log-change in price less than the reference value
-
- a = zeros( 1 , J )
- a[ I_Select ] = 1 # select those cases where the abs log-change in price is less than the 3/5 of the empirical cumulative density...
-
- A = a
- b = .5 # ... and assign the probability of these cases occuring as 50%. This moves the media of the distribution
-
+ Aeq <- matrix(1, 1, J)
+ beq <- 1
+
+ # ...constrain the median...
+ # absolute value of the log of changes in MSFT close prices
+ # (definition of realized volatility)
+ V <- abs(X[, 1])
+ # sorting of the abs value of log changes in prices from smallest to largest
+ V_Sort <- sort(V, decreasing = FALSE)
+ I_Sort <- order(V)
+
+ # represents the cumulative sum of probabilities from ~0 to 1
+ F <- cumsum(p[I_Sort])
+ # finds the (max) index corresponding to element with value <= 3/5 along the
+ # empirical cumulative density function for the abs log-changes in price
+ I_Reference <- max(matlab:::find(F <= 3 / 5))
+ # returns the corresponding abs log of change in price at the 3/5 of the
+ # cumulative density function
+ V_Reference <- V_Sort[I_Reference]
+
+ # finds all indices with value of abs log-change in price < reference value
+ I_Select <- find(V <= V_Reference)
+
+ a <- zeros(1, J)
+ # select those cases where the abs log-change in price is less than the 3/5 of
+ # the empirical cumulative density...
+ a[I_Select] <- 1
+
+ A <- a
+ # ... and assign the probability of these cases occuring as 50%.
+ # This moves the media of the distribution
+ b <- .5
+
# ...compute posterior probabilities
- p_ = EntropyProg( p , A , b , Aeq , beq )$p_
-
- return( p_ )
+ p_ <- EntropyProg(p, A, b, Aeq, beq)$p_
+
+ return(p_)
}
-#' Process views for the expectations and binding constraints as it appears in A. Meucci,
-#' "Fully Flexible Views: Theory and Practice", The Risk Magazine, October 2008,
-#' p 100-106
+#' @title Process views for the expectations and binding constraints
+#'
+#' @description Process views for the expectations and binding constraints as it
+#' appears in A. Meucci, "Fully Flexible Views: Theory and Practice", The Risk
+#' Magazine, October 2008, p 100-106.
#'
#' @param X : panel of joint factors realizations
#' @param p : vector of probabilities
@@ -242,42 +294,44 @@
#' @return p_ : vector of posterior probabilities
#'
#' @references
-#' A. Meucci, "Fully Flexible Views: Theory and Practice" \url{http://www.symmys.com/node/158}
+#' A. Meucci, "Fully Flexible Views: Theory and Practice"
+#' \url{http://www.symmys.com/node/158}
+#'
#' See Meucci script for "ButterflyTrading/ViewCurveSlope.m"
#'
-#' @author Xavier Valls \email{flamejat@@gmail.com}
+#' @author Xavier Valls \email{xaviervallspla@@gmail.com}
#' @export
-ViewCurveSlope = function( X , p )
-{
- # view 3
-
- J = nrow( X );
- K = ncol( X );
-
+ViewCurveSlope <- function(X, p) {
+ # view 3
+ J <- nrow(X)
+
# constrain probabilities to sum to one...
- Aeq = matrix( 1, 1 , J );
- beq = 1;
-
+ Aeq <- matrix(1, 1, J);
+ beq <- 1;
+
# ...constrain the expectation...
- V = X[ , 14 ] - X[ , 13 ];
- v = 0.0005;
-
- Aeq = rbind( Aeq , t(V) );
-
- beq = rbind( beq , v );
-
- A = b = matrix( nrow = 0 , ncol = 0 );
-
+ V <- X[, 14] - X[, 13];
+ v <- 0.0005;
+
+ Aeq <- rbind(Aeq, t(V));
+
+ beq <- rbind(beq, v);
+
+ A <- b <- matrix(nrow = 0, ncol = 0);
+
# ...compute posterior probabilities
- p_ = EntropyProg( p , A , b , Aeq ,beq )$p_;
+ p_ <- EntropyProg(p, A, b, Aeq,beq)$p_;
- return( p_ );
+ return(p_);
}
-#' Computes the conditional value at risk as it appears in A. Meucci, "Fully Flexible Views: Theory and Practice",
-#' The Risk Magazine, October 2008, p 100-106
-#'
+#' @title Computes the conditional value at risk
+#'
+#' @description Computes the conditional value at risk as it appears in
+#' A. Meucci, "Fully Flexible Views: Theory and Practice", The Risk Magazine,
+#' October 2008, p 100-106.
+#'
#' @param Units panel of joint factors realizations
#' @param Scenarios vector of probabilities
#' @param Conf Confidence
@@ -285,132 +339,150 @@
#' @return CVaR Conditional Value at Risk
#'
#' @references
-#' A. Meucci, "Fully Flexible Views: Theory and Practice" \url{http://www.symmys.com/node/158}
+#' A. Meucci, "Fully Flexible Views: Theory and Practice"
+#' \url{http://www.symmys.com/node/158}
+#'
#' See Meucci script for "ButterflyTrading/ComputeCVaR.m"
#'
#' @author Ram Ahluwalia \email{ram@@wingedfootcapital.com}
#' @export
-ComputeCVaR = function( Units , Scenarios , Conf )
+ComputeCVaR <- function(Units, Scenarios, Conf)
{
- PnL = Scenarios %*% Units
- Sort_PnL = PnL[ order( PnL , decreasing = FALSE ) ]
-
- J = length( PnL )
- Cut = round( J %*% ( 1 - Conf ) , 0 )
-
- CVaR = -mean( Sort_PnL[ 1:Cut ] )
-
- return( CVaR )
+ PnL <- Scenarios %*% Units
+ Sort_PnL <- PnL[order(PnL, decreasing = FALSE)]
+
+ J <- length(PnL)
+ Cut <- round(J %*% (1 - Conf), 0)
+
+ CVaR <- -mean(Sort_PnL[1:Cut])
+
+ return(CVaR)
}
-#' Computes the long-short conditional value at risk frontier as it appears in A. Meucci,
-#' "Fully Flexible Views: Theory and Practice", The Risk Magazine, October 2008, p 100-106
-#'
+#' @title Computes the long-short conditional value at risk
+#'
+#' @description Computes the long-short conditional value at risk frontier as it
+#' appears in A. Meucci, "Fully Flexible Views: Theory and Practice", The Risk
+#' Magazine, October 2008, p 100-106.
+#'
#' @param PnL Profit and Loss scenarios
#' @param Probs vector of probabilities
#' @param Butterflies list of securities with some analytics computed.
#' @param Options list of options
#'
#' @return Exp vector of expected returns for each asset
-#' @return SDev vector of security volatilities along the efficient frontier
+#' @return SDev vector of security volatilities along the efficient
+#' frontier
#' @return CVaR Conditional Value at Risk for each portfolio
-#' @return Composition matrix of compositions (security weights) for each portfolio along the efficient frontier
+#' @return Composition matrix of compositions (security weights) for each
+#' portfolio along the efficient frontier
#'
#' @references
-#' A. Meucci, "Fully Flexible Views: Theory and Practice" \url{http://www.symmys.com/node/158}
+#' A. Meucci, "Fully Flexible Views: Theory and Practice"
+#' \url{http://www.symmys.com/node/158}
+#'
#' See Meucci script for "ButterflyTrading/LongShortMeanCVaRFrontier.m"
#'
-#' @author Ram Ahluwalia \email{ram@@wingedfootcapital.com}, Xavier Valls \email{flamejat@@gmail.com}
+#' @author Ram Ahluwalia \email{ram@@wingedfootcapital.com},
+#' Xavier Valls \email{xaviervallspla@@gmail.com}
#' @export
-LongShortMeanCVaRFrontier = function( PnL , Probs , Butterflies , Options )
-{
- library( matlab )
- library( quadprog )
- library( limSolve )
-
+LongShortMeanCVaRFrontier <- function(PnL, Probs, Butterflies, Options) {
# setup constraints
- J = nrow(PnL); N = ncol(PnL)
- P_0s = matrix( , nrow = 1 , ncol = 0 )
- D_s = matrix( , nrow = 1 , ncol = 0 )
- emptyMatrix = matrix( nrow = 0 , ncol = 0 )
-
- for ( n in 1:N )
- {
- P_0s = cbind( P_0s , Butterflies[[n]]$P_0 ) # 1x9 matrix
- D_s = cbind( D_s , Butterflies[[n]]$Delta ) # 1x9 matrix
+ N <- ncol(PnL)
+ P_0s <- matrix(, nrow = 1, ncol = 0)
+ D_s <- matrix(, nrow = 1, ncol = 0)
+ emptyMatrix <- matrix(nrow = 0, ncol = 0)
+
+ for (n in 1:N) {
+ P_0s <- cbind(P_0s, Butterflies[[n]]$P_0) # 1x9 matrix
+ D_s <- cbind(D_s, Butterflies[[n]]$Delta) # 1x9 matrix
}
-
- Constr = list()
- Constr$Aeq = P_0s # linear coefficients in the constraints Aeq*X = beq (equality constraints)
- Constr$beq = Options$Budget # the constant vector in the constraints Aeq*x = beq
-
- if ( Options$DeltaNeutral == TRUE )
- {
- Constr$Aeq = rbind( Constr$Aeq , D_s ) # 2x9 matrix
- Constr$beq = rbind( Constr$beq , 0 ) # 2x9 matrix
+
+ Constr <- list()
+ # linear coefficients in the constraints Aeq*X = beq (equality constraints)
+ Constr$Aeq <- P_0s
+ # the constant vector in the constraints Aeq*x = beq
+ Constr$beq <- Options$Budget
+
+ if (Options$DeltaNeutral == TRUE) {
+ Constr$Aeq <- rbind(Constr$Aeq, D_s) # 2x9 matrix
+ Constr$beq <- rbind(Constr$beq, 0) # 2x9 matrix
}
-
- Constr$Aleq = rbind( diag( as.vector( P_0s ) ) , -diag( as.vector( P_0s ) ) ) # linear coefficients in the constraints A*x <= b. an 18x9 matrix
- Constr$bleq = rbind( Options$Limit * matrix( 1,N,1) , Options$Limit * matrix( 1,N,1) ) # constant vector in the constraints A*x <= b. an 18x1 matrix
-
+
+# linear coefficients in the constraints A*x <= b. an 18x9 matrix
+ Constr$Aleq <- rbind(diag(as.vector(P_0s)), -diag(as.vector(P_0s)))
+ # constant vector in the constraints A*x <= b. an 18x1 matrix
+ Constr$bleq <- rbind(Options$Limit * matrix(1,N,1),
+ Options$Limit * matrix(1,N,1))
+
# determine expectation of minimum-variance portfolio
- Exps = t(PnL) %*% Probs
- Scnd_Mom = t(PnL) %*% (PnL * (Probs %*% matrix( 1,1,N) ) )
- Scnd_Mom = ( Scnd_Mom + t(Scnd_Mom) ) / 2
- Covs = Scnd_Mom - Exps %*% t(Exps)
-
- Amat = rbind( Constr$Aeq , Constr$Aleq ) # stack the equality constraints on top of the inequality constraints
- bvec = rbind( Constr$beq , Constr$bleq ) # stack the equality constraints on top of the inequality constraints
-
- #if ( nrow(Covs) != length( zeros(N,1) ) ) stop("Dmat and dvec are incompatible!")
- #if ( nrow(Covs) != nrow(Amat)) stop("Amat and dvec are incompatible!")
-
- MinSDev_Units = solve.QP( Dmat = Covs , dvec = -1 * zeros(N,1) , Amat = -1*t(Amat) , bvec = -1*bvec , meq = length( Constr$beq) ) # TODO: Check this
- MinSDev_Exp = t( MinSDev_Units$solution ) %*% Exps
-
+ Exps <- t(PnL) %*% Probs
+ Scnd_Mom <- t(PnL) %*% (PnL * (Probs %*% matrix(1,1,N)))
+ Scnd_Mom <- (Scnd_Mom + t(Scnd_Mom)) / 2
+ Covs <- Scnd_Mom - Exps %*% t(Exps)
+
+# stack the equality constraints on top of the inequality constraints
+ Amat <- rbind(Constr$Aeq, Constr$Aleq)
+ bvec <- rbind(Constr$beq, Constr$bleq)
+
+ #if (nrow(Covs) != length(zeros(N,1))) stop("Dmat and dvec are incompatible!")
+ #if (nrow(Covs) != nrow(Amat)) stop("Amat and dvec are incompatible!")
+
+ MinSDev_Units <- solve.QP(Dmat = Covs, dvec = -1 * zeros(N,1),
+ Amat = -1 * t(Amat), bvec = -1 * bvec,
+ meq = length(Constr$beq)) # TODO: Check this
+ MinSDev_Exp <- t(MinSDev_Units$solution) %*% Exps
+
# determine expectation of maximum-expectation portfolio
-
- MaxExp_Units = linp( E = Constr$Aeq , F = Constr$beq , G = -1*Constr$Aleq , H = -1*Constr$bleq , Cost = -Exps , ispos = FALSE )$X
-
- MaxExp_Exp = t( MaxExp_Units ) %*% Exps
-
+
+ MaxExp_Units <- linp(E = Constr$Aeq, F = Constr$beq, G = -1 * Constr$Aleq,
+ H = -1 * Constr$bleq, Cost = -Exps, ispos = FALSE)$X
+
+ MaxExp_Exp <- t(MaxExp_Units) %*% Exps
+
# slice efficient frontier in NumPortf equally thick horizontal sections
- Grid = t( seq( from = Options$FrontierSpan[1] , to = Options$FrontierSpan[2] , length.out = Options$NumPortf ) )
- TargetExp = as.numeric( MinSDev_Exp ) + Grid * as.numeric( ( MaxExp_Exp - MinSDev_Exp ) )
-
+ Grid <- t(seq(from = Options$FrontierSpan[1], to = Options$FrontierSpan[2],
+ length.out = Options$NumPortf))
+ TargetExp <- as.numeric(MinSDev_Exp) + Grid * as.numeric(MaxExp_Exp -
+ MinSDev_Exp)
+
# compute composition, expectation, s.dev. and CVaR of the efficient frontier
- Composition = matrix( , ncol = N , nrow = 0 )
- Exp = matrix( , ncol = 1 , nrow = 0 )
- SDev = matrix( , ncol = 1 , nrow = 0 )
- CVaR = matrix( , ncol = 1 , nrow = 0 )
-
- for (i in 1:Options$NumPortf )
- {
- # determine least risky portfolio for given expectation
- AEq = rbind( Constr$Aeq , t(Exps) ) # equality constraint: set expected return for each asset...
- bEq = rbind( Constr$beq , TargetExp[i] )
-
- Amat = rbind( AEq , Constr$Aleq ) # stack the equality constraints on top of the inequality constraints
- bvec = rbind( bEq , Constr$bleq ) # ...and target portfolio return for i'th efficient portfolio
-
- # Why is FirstDegree "expected returns" set to 0?
+ Composition <- matrix(, ncol = N, nrow = 0)
+ Exp <- matrix(, ncol = 1, nrow = 0)
+ SDev <- matrix(, ncol = 1, nrow = 0)
+ CVaR <- matrix(, ncol = 1, nrow = 0)
+
+ for (i in 1:Options$NumPortf) {
+ ## determine least risky portfolio for given expectation
+ # equality constraint: set expected return for each asset...
+ AEq <- rbind(Constr$Aeq, t(Exps))
+ bEq <- rbind(Constr$beq, TargetExp[i])
+
+ # stack the equality constraints on top of the inequality constraints
+ Amat <- rbind(AEq, Constr$Aleq)
+ # ...and target portfolio return for i'th efficient portfolio
+ bvec <- rbind(bEq, Constr$bleq)
+
+ # Why is FirstDegree "expected returns" set to 0?
# Becasuse we capture the equality view in the equality constraints matrix
- # In other words, we have a constraint that the Expected Returns by Asset %*% Weights = Target Return
- Units = solve.QP( Dmat = Covs , dvec = -1*zeros(N,1) , Amat = -1*t(Amat) , bvec = -1*bvec , meq = length( bEq ) )
-
+ # In other words, we have a constraint that the Expected Returns by
+ # Asset %*% Weights = Target Return
+ Units <- solve.QP(Dmat = Covs, dvec = -1 * zeros(N,1), Amat = -1 * t(Amat),
+ bvec = -1 * bvec, meq = length(bEq))
+
# store results
- Composition = rbind( Composition , t( Units$solution ) )
-
- Exp = rbind( Exp , t( Units$solution ) %*% Exps )
- SDev = rbind( SDev , sqrt( t( Units$solution ) %*% Covs %*% Units$solution ) )
- CVaR = rbind( CVaR , ComputeCVaR( Units$solution , PnL , Options$Quant ) )
- }
-
- colnames( Composition ) = c( "MSFT_vol_30" , "MSFT_vol_91" , "MSFT_vol_182" ,
- "YHOO_vol_30" , "YHOO_vol_91" , "YHOO_vol_182" ,
- "GOOG_vol_30" , "GOOG_vol_91" , "GOOG_vol_182" )
-
- return( list( Exp = Exp , SDev = SDev , CVaR = CVaR , Composition = Composition ) )
+ Composition <- rbind(Composition, t(Units$solution))
+
+ Exp <- rbind(Exp, t(Units$solution) %*% Exps)
+ SDev <- rbind(SDev, sqrt(t(Units$solution) %*% Covs %*% Units$solution))
+ CVaR <- rbind(CVaR, ComputeCVaR(Units$solution, PnL, Options$Quant))
+ }
+
+ colnames(Composition) <- c("MSFT_vol_30", "MSFT_vol_91", "MSFT_vol_182",
+ "YHOO_vol_30", "YHOO_vol_91", "YHOO_vol_182",
+ "GOOG_vol_30", "GOOG_vol_91", "GOOG_vol_182")
+
+ return(list(Exp = Exp, SDev = SDev, CVaR = CVaR, Composition = Composition))
}
Modified: pkg/Meucci/R/EntropyProg.R
===================================================================
--- pkg/Meucci/R/EntropyProg.R 2015-06-16 14:22:14 UTC (rev 3676)
+++ pkg/Meucci/R/EntropyProg.R 2015-06-16 15:39:42 UTC (rev 3677)
@@ -98,11 +98,11 @@
if (((.999999 < sum(p)) & (sum(p) < 1.00001)) == FALSE)
stop("sum of probabilities from prior distribution must equal 1")
- if (nrow(Aeq) != nrow(beq)){
+ if (nrow(Aeq) != length(beq)){
stop("number of inequality constraints in matrix Aeq must match
number of elements in vector beq")
}
- if (nrow(A) != nrow(b)) {
+ if (nrow(A) != length(b)) {
stop("number of equality constraints in matrix A must match number of
elements in vector b")
}
Modified: pkg/Meucci/demo/ButterflyTrading.R
===================================================================
--- pkg/Meucci/demo/ButterflyTrading.R 2015-06-16 14:22:14 UTC (rev 3676)
+++ pkg/Meucci/demo/ButterflyTrading.R 2015-06-16 15:39:42 UTC (rev 3677)
@@ -1,101 +1,117 @@
-#' This script performs the butterfly-trading case study for the Entropy-Pooling approach by Attilio Meucci,
-#' as it appears in A. Meucci, "Fully Flexible Views: Theory and Practice", The Risk Magazine, October 2008,
-#' p 100-106
+#' This script performs the butterfly-trading case study for the Entropy-Pooling
+#'' approach by Attilio Meucci, as it appears in A. Meucci, "Fully Flexible
+#' Views: Theory and Practice", The Risk Magazine, October 2008, p 100-106
#'
#' Most recent version of article and MATLAB code available at
#' http://www.symmys.com/node/158
#'
#' @references
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/returnanalytics -r 3677
More information about the Returnanalytics-commits
mailing list