[Returnanalytics-commits] r3072 - in pkg/Meucci: . R data demo man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Sep 12 21:10:50 CEST 2013
Author: xavierv
Date: 2013-09-12 21:10:50 +0200 (Thu, 12 Sep 2013)
New Revision: 3072
Added:
pkg/Meucci/R/ButterflyTradingFunctions.R
pkg/Meucci/R/RankingInformationFunctions.R
pkg/Meucci/data/factorsDistribution.rda
pkg/Meucci/man/HorizonPricing.Rd
pkg/Meucci/man/ViewCurveSlope.Rd
pkg/Meucci/man/ViewImpliedVol.Rd
pkg/Meucci/man/ViewRealizedVol.Rd
pkg/Meucci/man/factorsDistribution.Rd
Modified:
pkg/Meucci/DESCRIPTION
pkg/Meucci/NAMESPACE
pkg/Meucci/R/BlackScholesCallPrice.R
pkg/Meucci/R/data.R
pkg/Meucci/data/butterfliesAnalytics.rda
pkg/Meucci/demo/ButterflyTrading.R
pkg/Meucci/man/BlackScholesCallPrice.Rd
pkg/Meucci/man/butterfliesAnalytics.Rd
pkg/Meucci/man/returnsDistribution.Rd
Log:
- documented some functions for the ButterflyAnalitics paper and changed the datafiles
Modified: pkg/Meucci/DESCRIPTION
===================================================================
--- pkg/Meucci/DESCRIPTION 2013-09-12 18:10:48 UTC (rev 3071)
+++ pkg/Meucci/DESCRIPTION 2013-09-12 19:10:50 UTC (rev 3072)
@@ -41,7 +41,6 @@
MASS,
reshape2,
Hmisc,
- fOptions,
moments,
nloptr,
ggplot2,
Modified: pkg/Meucci/NAMESPACE
===================================================================
--- pkg/Meucci/NAMESPACE 2013-09-12 18:10:48 UTC (rev 3071)
+++ pkg/Meucci/NAMESPACE 2013-09-12 19:10:50 UTC (rev 3072)
@@ -1,5 +1,7 @@
export(BlackLittermanFormula)
export(BlackScholesCallPrice)
+export(BlackScholesCallPutPrice)
+export(BlackScholesPutPrice)
export(Central2Raw)
export(CentralAndStandardizedStatistics)
export(CMAcombination)
@@ -25,6 +27,7 @@
export(GenerateLogNormalDistribution)
export(GenerateUniformDrawsOnUnitSphere)
export(hermitePolynomial)
+export(HorizonPricing)
export(integrateSubIntervals)
export(InterExtrapolate)
export(LeastInfoKernel)
@@ -64,3 +67,6 @@
export(SummStats)
export(Tweak)
export(TwoDimEllipsoid)
+export(ViewCurveSlope)
+export(ViewImpliedVol)
+export(ViewRealizedVol)
Modified: pkg/Meucci/R/BlackScholesCallPrice.R
===================================================================
--- pkg/Meucci/R/BlackScholesCallPrice.R 2013-09-12 18:10:48 UTC (rev 3071)
+++ pkg/Meucci/R/BlackScholesCallPrice.R 2013-09-12 19:10:50 UTC (rev 3072)
@@ -1,5 +1,5 @@
-#' Compute the Black-Scholes price of a European call option
-#' as described in A. Meucci, "Risk and Asset Allocation", Springer, 2005.
+#' Compute the Black-Scholes price of a European call or put option
+#' as described in A. Meucci, "Risk and Asset Allocation", Springer, 2005.
#'
#' @param spot : [scalar] spot price of underlying
#' @param K : [scalar] strike of the call optioon
@@ -8,7 +8,8 @@
#' @param T : [scalar] time to maturity in years
#'
#' @return c : [scalar] price of European call(s)
-#' @return delta : [scalar] delta of the call(s)
+#' @return p : [scalar] price of European put(s)
+#' @return delta : [scalar] delta of the call(s) or put(s)
#' @return cash : [scalar] cash held in a replicating portfolio
#'
#' @note
@@ -21,13 +22,41 @@
#' @author Xavier Valls \email{flamejat@@gmail.com}
#' @export
-BlackScholesCallPrice = function(spot, K, r, vol, T)
+BlackScholesCallPrice = function( spot, K, r, vol, T )
{
- d1 = ( log( spot / K ) + ( r + vol * vol / 2) * T) / (vol * sqrt(T));
- d2 = d1 - vol * sqrt(T);
+ d1 = ( log( spot / K ) + ( r + vol * vol / 2) * T) / (vol * sqrt(T));
+ d2 = d1 - vol * sqrt(T);
delta = pnorm(d1);
- cash = -K * exp( -r * T ) * pnorm( d2 );
- c = spot * delta + cash;
+ cash = -K * exp( -r * T ) * pnorm( d2 );
+ c = spot * delta + cash;
return( list( c = c, delta = delta, cash = cash ) );
-}
\ No newline at end of file
+}
+
+#' @rdname BlackScholesCallPrice
+#' @export
+
+BlackScholesPutPrice = function( spot, K, r, vol, T )
+{
+ d1 = ( log( spot / K ) + ( r + vol * vol / 2) * T) / (vol * sqrt(T));
+ d2 = d1 - vol * sqrt(T);
+ delta = pnorm( -d1 );
+ cash = -K * exp( -r * T ) * pnorm( d2 );
+ p = -( spot * delta + cash );
+
+ return( list( put = p, delta = delta, cash = cash ) );
+}
+
+#' @rdname BlackScholesCallPrice
+#' @export
+
+BlackScholesCallPutPrice = function( spot, K, r, vol, T )
+{
+ d1 = ( log( spot / K ) + ( r + vol * vol / 2) * T) / (vol * sqrt(T));
+ d2 = d1 - vol * sqrt(T);
+ cash = -K * exp( -r * T ) * pnorm( d2 );
+ c = spot * pnorm( d1 ) + cash;
+ p = -( spot * pnorm( -d1 ) + cash);
+
+ return( list( call = c, put = p, cash = cash ) );
+}
Added: pkg/Meucci/R/ButterflyTradingFunctions.R
===================================================================
--- pkg/Meucci/R/ButterflyTradingFunctions.R (rev 0)
+++ pkg/Meucci/R/ButterflyTradingFunctions.R 2013-09-12 19:10:50 UTC (rev 3072)
@@ -0,0 +1,394 @@
+# 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 )
+}
+
+#' 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
+#'
+#' @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"
+#'
+#' @author Ram Ahluwalia \email{ram@@wingedfootcapital.com} and Xavier Valls \email{flamejat@@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'
+
+ # 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 ] }
+
+ # 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 #################################
+ #
+ ## 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 )
+ #
+ ####################################################################################################
+
+ 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 )
+}
+
+ViewCurveSlopeTest = function( X , p )
+{
+ J = nrow( X ) ; K = ncol( X )
+
+ # constrain probabilities to sum to one...
+ Aeq = matrix( 1, 1 , J )
+ beq = matrix( 1 , nrow = 1 , ncol = 1 )
+ browser()
+ # ...constrain the expectation...
+ 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 )
+
+ # 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
+
+ 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
+
+ # ...compute posterior probabilities
+ p_ = EntropyProg( p , A , b , Aeq ,beq )
+ 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.
+#'
+#' @param X : panel of joint factors realizations
+#' @param p : vector of probabilities
+#'
+#' @return p_ : vector of posterior probabilities
+#'
+#' @references
+#' 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}
+#' @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 );
+
+ # constrain probabilities to sum to one...
+ Aeq = matrix( 1, 1 , J )
+ beq = 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
+
+ # ...compute posterior probabilities
+ 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
+#'
+#' @param X : panel of joint factors realizations
+#' @param p : vector of probabilities
+#'
+#' @return p_ : vector of posterior probabilities
+#'
+#' @references
+#' 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}
+#' @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 )
+
+ # 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
+
+ # ...compute posterior probabilities
+ 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
+#'
+#' @param X : panel of joint factors realizations
+#' @param p : vector of probabilities
+#'
+#' @return p_ : vector of posterior probabilities
+#'
+#' @references
+#' 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}
+#' @export
+
+ViewCurveSlope = function( X , p )
+{
+ # view 3
+
+ J = nrow( X );
+ K = ncol( X );
+
+ # constrain probabilities to sum to one...
+ 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 );
+
+ # ...compute posterior probabilities
+ p_ = EntropyProg( p , A , b , Aeq ,beq )$p_;
+
+ return( p_ );
+}
+
+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 )
+}
+
+LongShortMeanCVaRFrontier = function( PnL , Probs , Butterflies , Options )
+{
+ library( matlab )
+ library( quadprog )
+ library( limSolve )
+
+ # 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
+ }
+
+ 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$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
+
+ # 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
+
+ # 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
+
+ # 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 ) )
+
+ # 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?
+ # 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 ) )
+
+ # 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 ) )
+}
+
+
+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 )
+}
+
Added: pkg/Meucci/R/RankingInformationFunctions.R
===================================================================
--- pkg/Meucci/R/RankingInformationFunctions.R (rev 0)
+++ pkg/Meucci/R/RankingInformationFunctions.R 2013-09-12 19:10:50 UTC (rev 3072)
@@ -0,0 +1,224 @@
+# TODO: add max weights constraint to EfficientFrontier()
+# TODO: add computeCVaR to EfficientFrontier()
+# TODO: confirm QuadProg does not have a bug (i.e. it can optimize expected returns without use dvec by adding an equality constraint)
+
+#' Plots the efficient frontier, as it appears in A. Meucci, "Fully Flexible Views: Theory and Practice", The Risk Magazine,
+#' October 2008, p 100-106.
+#'
+#' @param e the NumPortf x 1 matrix of expected returns for each portfolio along the efficient frontier
+#' @param s the NumPortf x 1 matrix of standard deviation of returns for each portfolio along the efficient frontier
+#' @param w the NumPortf x N 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}
+#' See Meucci script for "RankingInformation/PlotFrontier.m"
+#'
+#' @author Xavier Valls \email{flamejat@@gmail.com}
+#' @export
+
+
+PlotFrontier = function( e, s, w )
+{
+ xx = dim( w )[ 1 ];
+ N = dim( w )[ 2 ];
+ Data = t( apply( w, 1, cumsum ) );
+
+ plot( c(min(s), 0), xlim = c( min(s) , max(s) ), ylim = c( 0, max(Data) ),
+ main= "frontier", xlab = " Portfolio # risk propensity", ylab = "Portfolio composition" );
+
+ for( n in 1 : N )
+ {
+ x = rbind( min(s), s, max(s) );
+ y = rbind( 0, matrix( Data[ , N-n+1 ] ), 0 );
+ polygon( x, y, col = rgb( 0.9 - mod(n,3)*0.2, 0.9 - mod(n,3)*0.2, 0.9 - mod(n,3)*0.2) );
+ }
+}
+
+#' Plots the results of computing the efficient frontier (Expected returns and frontier), as it appears in A. Meucci, "Fully Flexible Views: Theory and Practice", The Risk Magazine,
+#' October 2008, p 100-106.
+#'
+#' @param e the NumPortf x 1 matrix of expected returns for each portfolio along the efficient frontier
+#' @param s the NumPortf x 1 matrix of standard deviation of returns for each portfolio along the efficient frontier
+#' @param w the NumPortf x N matrix of compositions (security weights) for each portfolio along the efficient frontier
+#' @param M the NumPortf x 1 vector of expected returns for each asset
+#' @param Lower constraints
+#' @param Upper constraints
+#'
+#' @references
+#' A. Meucci, "Fully Flexible Views: Theory and Practice" \url{http://www.symmys.com/node/158}
+#' See Meucci script for "RankingInformation/PlotResults.m"
+#'
+#' @author Xavier Valls \email{flamejat@@gmail.com}
+
+PlotResults = function( e, s, w, M, Lower = NULL , Upper = NULL)
+{
+ N = length( M );
+ dev.new();
+ par( mfrow = c( 1, 2 ) );
+ h1 = hist( M*100, plot = F )
+ barplot( h1$density, horiz = T, main = "expected returns", xlab = "", ylab = "" );
+ if(length(Lower) || length(Upper))
+ {
+ Changed = array( 0, N );
+ Changed[ union( Lower, Upper ) ] = M[ union( Lower, Upper ) ] * 100;
+ h2 = hist(Changed, plot = F );
+ barplot( h2$density, horiz = T, col = "red", add = T );
+ }
+
+ PlotFrontier( e*100, s*100, w );
+}
+
+
+
+#' Computes posterior probabilities to view the rankings, as it appears in A. Meucci,
+#' "Fully Flexible Views: Theory and Practice", The Risk Magazine, October 2008, p 100-106.
+#'
+#' @param X a vector containing returns for all the asset classes
+#' @param p a vector containing the prior probability values
+#' @param Lower a vector of indexes indicating which column is lower than the corresponding column number in Upper
+#' @param Upper a vector of indexes indicating which column is lower than the corresponding column number in Upper
+#'
+#' @references
+#' A. Meucci, "Fully Flexible Views: Theory and Practice" \url{http://www.symmys.com/node/158}
+#' See Meucci script for "RankingInformation/ViewRanking.m"
+#'
+#' @author Ram Ahluwalia \email{ram@@wingedfootcapital.com}
+#' @export EntropyProg
+
+# example ViewRanking( X , p , Lower = c(3,4) , Upper = c(4,5) ) # two inequality views: asset 3 < asset 4 returns, and asset 4 < asset 5 returns
+
+ViewRanking = function( X , p , Lower , Upper )
+{
+ library( matlab )
+ J = nrow( X )
+ N = ncol( X )
+
+ K = length( Lower )
+
+ # constrain probabilities to sum to one across all scenarios...
+ Aeq = ones( 1 , J )
+ beq = 1
+
+ # ...constrain the expectations... A*x <= 0
+ # X[,Lower] refers to the column of returns for Asset-lower
+ # X[,Upper] refers to the column of returns for Asset-lower
+ # X[ , Lower ] - X[ , Upper ] is vector returns of the "lower"" asset less the returns of the "higher" asset
+ V = X[ , Lower ] - X[ , Upper ] # Jx1 vector. Expectation is assigned to each scenario
+
+ A = t( V )
+ b = 0 # The expectation is that (Lower - Upper)x <= 0. (i.e. The returns of upper are greater than zero for each scenario)
+
+ # ...compute posterior probabilities
+ p_ = EntropyProg( p , A , as.matrix(b) , Aeq , as.matrix(beq) )
+
+ return( p_ )
+}
+
+#' Generates an efficient frontier based on Meucci's Ranking Information version and returns a A list with
+#' NumPortf efficient portfolios whos returns are equally spaced along the whole range of the efficient frontier,
+#' 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
+#'
+#' @param X a matrix with the joint-scenario probabilities by asset (rows are joint-scenarios, columns are assets)
+#' @param p a vector of probabilities associated with each scenario in matrix X
+#' @param Options a list of options....TBD
+#'
+#' @return Exps the NumPortf x 1 vector of expected returns for each asset
+#' @return Covs the NumPortf x N vector of security volatilities along the efficient frontier
+#' @return w the NumPortf x N matrix of compositions (security weights) for each portfolio along the efficient frontier
+#' @return e the NumPortf x 1 matrix of expected returns for each portfolio along the efficient frontier
+#' @return s the NumPortf x 1 matrix of standard deviation of returns for each portfolio along the efficient frontier
+#'
+#' @references
+#' A. Meucci, "Fully Flexible Views: Theory and Practice" \url{http://www.symmys.com/node/158}
+#' See Meucci script for "RankingInformation/EfficientFrontier.m"
+#'
+#' @author Ram Ahluwalia \email{ram@@wingedfootcapital.com} and Xavier Valls \email{flamejat@@gmail.com}
+#' @export
+
+RIEfficientFrontier = function( X , p , Options)
+{
+
+ if( !require("limSolve") ) stop("This script requieres the limSolve package installed")
+
+
+ library( matlab )
+
+ J = nrow( X ) # number of scenarios
+ N = ncol( X ) # number of assets
+
+ Exps = t(X) %*% p # probability-weighted expected return of each asset
+
+ Scnd_Mom = t(X) %*% (X * ( p %*% matrix( 1, 1 , N ) ) )
+ Scnd_Mom = ( Scnd_Mom + t(Scnd_Mom) ) / 2 # an N*N matrix
+ Covs = Scnd_Mom - Exps %*% t( Exps )
+
+ Constr = list()
+
+ # constrain the sum of weights to 1
+ Constr$Aeq = matrix( 1, 1 , N )
+ Constr$beq = 1
+
+ # constrain the weight of any security to between 0 and 1
+ Constr$Aleq = rbind( diag( 1, N ) , - diag( 1, N ) ) # linear coefficients matrix A in the inequality constraint A*x <= b
+ Constr$bleq = rbind( matrix( 1, N, 1 ) , matrix( 0, N, 1 ) ) # constraint vector b in the inequality constraint A*x <= b
+
+ 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
+
+ ############################################################################################
+ # determine return of minimum-risk portfolio
+ FirstDegree = matrix( 0, N , 1 ) # TODO: assumes that securities have zero expected returns when computing efficient frontier?
+ SecondDegree = Covs
+ # Why is FirstDegree "expected returns" set to 0?
+ # 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
+ MinVol_Weights = solve.QP( Dmat = SecondDegree , dvec = -1*FirstDegree , Amat = -1*t(Amat) , bvec = -1*bvec , meq = length( Constr$beq ) )
+ MinSDev_Exp = t( MinVol_Weights$solution ) %*% Exps
+
+ ############################################################################################
+ # determine return of maximum-return portfolio
+ FirstDegree = -Exps
+ MaxRet_Weights = linp( E = Constr$Aeq , F = Constr$beq , G = -1*Constr$Aleq , H = -1*Constr$bleq , Cost = FirstDegree , ispos = FALSE )$X
+ MaxExp_Exp = t( MaxRet_Weights) %*% Exps
+
+ ############################################################################################
+ # slice efficient frontier in NumPortf equally thick horizontal sections
+ Grid = matrix( , ncol = 0 , nrow = 0 )
+ Grid = t( seq( from = Options$FrontierSpan[1] , to = Options$FrontierSpan[2] , length.out = Options$NumPortf ) )
+
+ # the portfolio return varies from a minimum of MinSDev_Exp up to a maximum of MaxExp_Exp
+ # We establish equally-spaced portfolio return targets and use this find efficient portfolios
+ # in the next step
+ Targets = as.numeric( MinSDev_Exp ) + Grid * as.numeric( ( MaxExp_Exp - MinSDev_Exp ) )
+
+ ############################################################################################
+ # compute the NumPortf compositions and risk-return coordinates
+ FirstDegree = matrix( 0, N , 1 )
+
+ w = matrix( , ncol = N , nrow = 0 )
+ e = matrix( , ncol = 1 , nrow = 0 )
+ s = matrix( , ncol = 1 , nrow = 0 )
+
+ for ( i in 1:Options$NumPortf )
+ {
+ # determine least risky portfolio for given expected return
+ # Ax = b ; Exps %*% weights = Target Return
+ AEq = rbind( Constr$Aeq , t( Exps ) ) # equality constraint: set expected return for each asset...
+ bEq = rbind( Constr$beq , Targets[ i ] ) # ...and target portfolio return for i'th efficient portfolio
+
+ Amat = rbind( AEq , Constr$Aleq ) # stack the equality constraints on top of the inequality constraints
+ bvec = rbind( bEq , Constr$bleq )
+
+ Weights = solve.QP( Dmat = SecondDegree , dvec = -1*FirstDegree , Amat = -1*t(Amat) , bvec = -1*bvec , meq = length( bEq ) )
+
+ w = rbind( w , Weights$solution )
+ s = rbind( s , sqrt( t(Weights$solution) %*% Covs %*% Weights$solution ) )
+ e = rbind( e , Weights$solution %*% Exps )
+ }
+
+ return( list( e = e , Sdev = s , Composition = w , Exps = Exps , Covs = Covs ) )
+}
Modified: pkg/Meucci/R/data.R
===================================================================
--- pkg/Meucci/R/data.R 2013-09-12 18:10:48 UTC (rev 3071)
+++ pkg/Meucci/R/data.R 2013-09-12 19:10:50 UTC (rev 3072)
@@ -210,7 +210,7 @@
#' @keywords data
NULL
-#' @title Panel X of joint returns realizations and vector p of respective probabilities
+#' @title Panel X of joint returns realizations and vector p of respective probabilities for returns
#'
#' @name returnsDistribution
#' @docType data
@@ -220,9 +220,9 @@
#' @keywords data
NULL
-#' @title Factor Distribution Butterflies
+#' @title Panel X of joint factors realizations and vector p of respective probabilities for factors
#'
-#' @name FDButterflies
+#' @name factorsDistribution
#' @docType data
#' @author Xavier Valls\email{flamejat@@gmail.com}
#' @references A. Meucci, "Fully Flexible Views: Theory and Practice", The Risk Magazine,
@@ -230,7 +230,7 @@
#' @keywords data
NULL
-#' @title Butterflies Analytics
+#' @title list of securities with analytics computed.
#'
#' @name butterfliesAnalytics
#' @docType data
Modified: pkg/Meucci/data/butterfliesAnalytics.rda
===================================================================
(Binary files differ)
Added: pkg/Meucci/data/factorsDistribution.rda
===================================================================
(Binary files differ)
Property changes on: pkg/Meucci/data/factorsDistribution.rda
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Modified: pkg/Meucci/demo/ButterflyTrading.R
===================================================================
--- pkg/Meucci/demo/ButterflyTrading.R 2013-09-12 18:10:48 UTC (rev 3071)
+++ pkg/Meucci/demo/ButterflyTrading.R 2013-09-12 19:10:50 UTC (rev 3072)
@@ -9,27 +9,25 @@
#' A. Meucci, "Fully Flexible Views: Theory and Practice" \url{http://www.symmys.com/node/158}
#' See Meucci script for "ButterflyTrading/S_MAIN.m"
#'
-#' @author Ram Ahluwalia \email{ram@@wingedfootcapital.com} and Xavier Valls \email{flamejat@@gmail.com}
+#' @author Xavier Valls \email{flamejat@@gmail.com} and Ram Ahluwalia \email{ram@@wingedfootcapital.com}
+###########################################################################################################
+# Load panel X of joint factors realizations and vector p of respective probabilities
+# In real life, these are provided by the estimation process
+###########################################################################################################
-###################################################################
-#' Load panel X of joint factors realizations and vector p of respective probabilities
-#' In real life, these are provided by the estimation process
-###################################################################
-load("butterflyTradingX.rda")
+load( "../data/factorsDistribution.rda" )
-#library( R.matlab )
-#library( matlab )
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/returnanalytics -r 3072
More information about the Returnanalytics-commits
mailing list