[Returnanalytics-commits] r1907 - in pkg/PerformanceAnalytics/sandbox: . Meucci
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Apr 23 22:30:34 CEST 2012
Author: braverock
Date: 2012-04-23 22:30:33 +0200 (Mon, 23 Apr 2012)
New Revision: 1907
Added:
pkg/PerformanceAnalytics/sandbox/Meucci/
pkg/PerformanceAnalytics/sandbox/Meucci/.~lock.Meucci_functions.csv#
pkg/PerformanceAnalytics/sandbox/Meucci/FactorDistributions.rda
pkg/PerformanceAnalytics/sandbox/Meucci/MeucciAnalyticalvsNumerical.R
pkg/PerformanceAnalytics/sandbox/Meucci/MeucciButterflyTrading.R
pkg/PerformanceAnalytics/sandbox/Meucci/MeucciCmaCopula.R
pkg/PerformanceAnalytics/sandbox/Meucci/MeucciEntropyProg.R
pkg/PerformanceAnalytics/sandbox/Meucci/MeucciFreaqEst.rda
pkg/PerformanceAnalytics/sandbox/Meucci/MeucciFullFlexibleBayesNets.R
pkg/PerformanceAnalytics/sandbox/Meucci/MeucciHermiteGrid.R
pkg/PerformanceAnalytics/sandbox/Meucci/MeucciInvariantProjection.R
pkg/PerformanceAnalytics/sandbox/Meucci/MeucciRankingInformation.R
pkg/PerformanceAnalytics/sandbox/Meucci/MeucciReturnsDistribution.rda
pkg/PerformanceAnalytics/sandbox/Meucci/MeucciRobustBayesianAllocation.R
pkg/PerformanceAnalytics/sandbox/Meucci/MeucciTweakTest.rda
pkg/PerformanceAnalytics/sandbox/Meucci/Meucci_DetectOutliersviaMVE.R
pkg/PerformanceAnalytics/sandbox/Meucci/Meucci_functions.csv
pkg/PerformanceAnalytics/sandbox/Meucci/butterflyTradingX.rda
pkg/PerformanceAnalytics/sandbox/Meucci/logToArithmeticCovariance.R
Log:
- add Meucci functions and data ported from original Matlab code contributed by Ram Ahluwalia
Added: pkg/PerformanceAnalytics/sandbox/Meucci/.~lock.Meucci_functions.csv#
===================================================================
--- pkg/PerformanceAnalytics/sandbox/Meucci/.~lock.Meucci_functions.csv# (rev 0)
+++ pkg/PerformanceAnalytics/sandbox/Meucci/.~lock.Meucci_functions.csv# 2012-04-23 20:30:33 UTC (rev 1907)
@@ -0,0 +1 @@
+Brian Peterson,brian,brian-rcg,23.04.2012 15:26,file:///home/brian/.libreoffice/3;
\ No newline at end of file
Added: pkg/PerformanceAnalytics/sandbox/Meucci/FactorDistributions.rda
===================================================================
(Binary files differ)
Property changes on: pkg/PerformanceAnalytics/sandbox/Meucci/FactorDistributions.rda
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: pkg/PerformanceAnalytics/sandbox/Meucci/MeucciAnalyticalvsNumerical.R
===================================================================
--- pkg/PerformanceAnalytics/sandbox/Meucci/MeucciAnalyticalvsNumerical.R (rev 0)
+++ pkg/PerformanceAnalytics/sandbox/Meucci/MeucciAnalyticalvsNumerical.R 2012-04-23 20:30:33 UTC (rev 1907)
@@ -0,0 +1,140 @@
+# TODO: Determine how to extend correlation view to multiple assets
+# TODO: Create plot distributions function
+
+PlotDistributions = function( X , p , Mu , Sigma , p_ , Mu_ , Sigma_ )
+ {
+ J = nrow( X )
+ N = ncol( X )
+
+ NBins = round( 10*log( J ) )
+
+ for ( n in 1:N )
+ {
+ # set ranges
+ xl = min(X[ , n ] )
+ xh = max(X[ , n ] )
+ #x = [xl : (xh-xl)/100 : xh]
+
+ # posterior numerical
+ h3 = pHist(X[ ,n] , p_ , NBins )
+
+ # posterior analytical
+ h4 = plot( x , normpdf( x , Mu_[n] , sqrt( Sigma_[n,n] ) ) )
+
+ # prior analytical
+ h2 = plot( x , normpdf( x , Mu[n] ,sqrt( Sigma[n,n] ) ) )
+
+ # xlim( cbind( xl , xh ) )
+ # legend([h3 h4 h2],'numerical', 'analytical', 'prior')
+ }
+ }
+
+
+#' Calculate the full-confidence posterior distributions of Mu and Sigma
+#'
+#' @param M a numeric vector with the Mu of the normal reference model
+#' @param Q a numeric vector used to construct a view on expectation of the linear combination Q %*% X
+#' @param M_Q a numeric vector with the view of the expectations of QX
+#' @param S a covariance matrix for the normal reference model
+#' @param G a numeric vector used to construct a view on covariance of the linear combination G %*% X
+#' @param S_G a numeric with the expectation associated with the covariance of the linear combination GX
+#'
+#' @return a list with
+#' M_ a numeric vector with the full-confidence posterior distribution of Mu
+#' S_ a covariance matrix with the full-confidence posterior distribution of Sigma
+#'
+#' @references
+#' \url{http://www.symmys.com}
+#' See Meucci script Prior2Posterior.m attached to Entropy Pooling Paper
+#' @author Ram Ahluwalia \email{ram@@wingedfootcapital.com}
+Prior2Posterior = function( M , Q , M_Q , S , G , S_G )
+ {
+ # See Appendix A.1 formula 49 for derivation
+ M_ = M + S %*% t(Q) %*% solve( Q %*% S %*% t(Q) ) %*% (M_Q - Q %*% M)
+
+ # See Appendix A.1 formula 57 for derivation
+ S_= S + (S %*% t(G)) %*% ( solve( G %*% S %*% t(G) ) %*% S_G %*% solve( G %*% S %*% t(G) ) - solve( G%*%S%*%t(G) ) ) %*% ( G %*% S )
+
+ return ( list( M_ = M_ , S_ = S_ ) )
+ }
+
+
+# This example script compares the numerical and the analytical solution of entropy-pooling, see
+# "A. Meucci - Fully Flexible Views: Theory and Practice
+# and associated script S_Main
+# Example compares analytical vs. numerical approach to entropy pooling
+
+# Code by A. Meucci, September 2008
+# Last version available at www.symmys.com > Teaching > MATLAB
+
+###############################################################
+# prior
+###############################################################
+library( matlab )
+library( MASS )
+# analytical representation
+N = 2 # market dimension (2 assets)
+Mu = zeros( N , 1 )
+r= .6
+Sigma = ( 1 - r ) * eye( N ) + r * ones( N , N ) # nxn correlation matrix with correlaiton 'r' in off-diagonals
+
+# numerical representation
+J = 100000 # number of scenarios
+p = ones( J , 1 ) / J
+dd = mvrnorm( J / 2 , zeros( N , 1 ) , Sigma ) # distribution centered on (0,0) with variance Sigma
+X = ones( J , 1 ) %*% t(Mu) + rbind( dd , -dd ) # JxN matrix of scenarios
+
+###############################################################
+# views
+###############################################################
+
+# location
+Q = matrix( c( 1 , -1 ) , nrow = 1 ) # long the first and asset and short the second asset produces an expectation (of Mu_Q calculated below)
+Mu_Q = .5
+
+# scatter
+G = matrix( c( -1 , 1 ) , nrow = 1 )
+Sigma_G = .5^2
+
+###############################################################
+# posterior
+###############################################################
+
+# analytical posterior
+RevisedMuSigma = Prior2Posterior( Mu , Q , Mu_Q , Sigma , G , Sigma_G )
+Mu_ = RevisedMuSigma$M_
+
+# numerical posterior
+Aeq = ones( 1 , J ) # constrain probabilities to sum to one...
+beq = 1
+
+# create views
+ QX = X %*% t(Q) # a Jx1 matrix
+
+ Aeq = rbind( Aeq , t(QX) ) # ...constrain the first moments...
+ # QX is a linear combination of vector Q and the scenarios X
+
+ beq = rbind( beq , Mu_Q )
+
+ SecMom = G %*% Mu_ %*% t(Mu_) %*% t(G) + Sigma_G # ...constrain the second moments...
+ # We use Mu_ from analytical result. We do not use Revised Sigma because we are testing whether
+ # the numerical approach for handling expectations of covariance matches the analytical approach
+ # TODO: Can we perform this procedure without relying on Mu_ from the analytical result?
+ GX = X %*% t(G)
+
+ for ( k in 1:nrow( G ) )
+ {
+ for ( l in k:nrow( G ) )
+ {
+ Aeq = rbind( Aeq , t(GX[ , k ] * GX[ , l ] ) )
+ beq = rbind( beq , SecMom[ k , l ] )
+ }
+ }
+
+emptyMatrix = matrix( , nrow = 0 , ncol = 0 )
+p_ = EntropyProg( p , emptyMatrix , emptyMatrix , Aeq , beq ) # ...compute posterior probabilities
+
+###############################################################
+# plots
+###############################################################
+PlotDistributions( X , p , Mu , Sigma , p_ , Mu_ , Sigma_ )
\ No newline at end of file
Added: pkg/PerformanceAnalytics/sandbox/Meucci/MeucciButterflyTrading.R
===================================================================
--- pkg/PerformanceAnalytics/sandbox/Meucci/MeucciButterflyTrading.R (rev 0)
+++ pkg/PerformanceAnalytics/sandbox/Meucci/MeucciButterflyTrading.R 2012-04-23 20:30:33 UTC (rev 1907)
@@ -0,0 +1,445 @@
+# Todo: plot efficient frontier in R
+
+PlotFrontier = function( e , s , w )
+ {
+ # subplot(2,1,1)
+ plot( s , e )
+ # grid on
+ # set(gca,'xlim',[min(s) max(s)])
+ #
+ # subplot(2,1,2)
+
+ xx = nrow( w ) ; N = ncol( w )
+ Data = apply( w , 1 , cumsum ) #TODO: Check. Take cumulative sum of *rows*. Try sapply?
+
+ for ( n in 1:N )
+ {
+ x = cbind( min(s) , s , max(s) )
+ y = cbind( 0 , Data[ , N-n+1 ] , 0 )
+ # hold on
+ #h = fill( x , y , cbind( .9 , .9 , .9) - mod( n , 3 ) %*% cbind( .2 , .2 , .2) )
+ }
+
+ #set(gca,'xlim',[min(s) max(s)],'ylim',[0 max(max(Data))])
+ #xlabel('portfolio # (risk propensity)')
+ #ylabel('portfolio composition')
+ }
+
+ViewCurveSlope = function( X , p )
+ {
+ # view 3 (expectations and binding constraints): slope of the yield curve will increase by 5 bp
+
+ J = nrow( X ) ; K = ncol( X )
+
+ # constrain probabilities to sum to one...
+ Aeq = ones( 1 , J )
+ beq = 1
+
+ # ...constrain the expectation...
+ V = X[ , 14 ] - X[ , 13 ]
+ v = .0005
+
+ Aeq = rbind( Aeq , t(V) )
+
+ beq = rbind( beq , v )
+
+ A = b = emptyMatrix
+
+ # ...compute posterior probabilities
+ p_ = EntropyProg( p , A , b , Aeq ,beq )$p_
+ return( p_ )
+ }
+
+ViewRealizedVol = function( X , p )
+ {
+ # view 2 (relative inequality view on median): 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 = ones( 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_ )
+ }
+
+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 = ones( 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_ )
+ }
+
+ComputeCVaR = function( Units , Scenarios , Conf )
+ {
+ PnL = Scenarios %*% Units
+ Sort_PnL = PnL[ order( PnL , decreasing = FALSE ) ] # DOUBLE CHECK IF I SHOULD USE ORDER INSTEAD OF SORT
+
+ 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 * ones(N,1) , Options$Limit * ones(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 %*% ones(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=-.00000000001
+ b= .00000000001
+
+ s = sig + a/sqrt(T) * ( log(K) - log(y) ) + b/T*( log(K) - log(y) )^2
+
+ return( s )
+ }
+
+HorizonPricing = function( Butterflies , X )
+{
+r = .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 )
+
+ # 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 )
+
+ Butterflies[[s]]$P_T = C at price + P at price
+ PnL = cbind( PnL , Butterflies[[s]]$P_T )
+ }
+ PnL = PnL[ , -1 ]
+}
+
+#' 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"
+#' available at www.symmys.com > Research > Working Papers
+#' Adapted from Code by A. Meucci, September 2008
+#' Last version available at www.symmys.com > Teaching > MATLAB
+
+
+###################################################################
+#' 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")
+
+FactorNames = c('MSFT_close' , 'MSFT_vol_30' , 'MSFT_vol_91' , 'MSFT_vol_182' ,
+ 'YHOO_close' , 'YHOO_vol_30' , 'YHOO_vol_91' , 'YHOO_vol_182' ,
+ 'GOOG_close' , 'GOOG_vol_30' , 'GOOG_vol_91' , 'GOOG_vol_182' ,
+ 'USD SWAP 2Y rate' , 'USD SWAP 10Y rate' )
+
+p = matrix( rep( 1 / 100224 , 100224 ) , ncol = 1 ) # creates a Jx1 matrix of equal probabilities for each scenario
+
+ans = 1 / 100224
+
+library( R.matlab )
+library( matlab )
+
+emptyMatrix = matrix( nrow = 0 , ncol = 0 )
+
+###################################################################
+# Load current prices, deltas and other analytics of the securities
+# In real life, these are provided by data provider
+###################################################################
+
+load( "FactorDistributions.rda" )
+
+# create Butterflies as a list of named arrays
+Butterflies = as.matrix( Butterflies[[1]] , nrow = 8 , ncol = 9 )
+Butterflies = matrix(Butterflies, ncol = 9 , nrow = 8 )
+rownames( Butterflies ) = c( "Name" , "P_0" , "Y_0" , "K" , "T" , "sig_0" , "Delta" , "Vega" )
+colnames( Butterflies ) = 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" )
+
+colnames( X ) = FactorNames
+Butterflies = lapply( seq_len( ncol( Butterflies ) ), function( i ) Butterflies[ , i ] )
+
+###################################################################
+# Map factors scenarios into p&l scenarios at the investment horizon
+# In real life with complex products, the pricing can be VERY costly
+###################################################################
+PnL = HorizonPricing( Butterflies , X )
+
+###################################################################
+# compute mean-risk efficient frontier based on prior market distribution
+###################################################################
+Options = list()
+Options$Quant = .95 # quantile level for CVaR
+Options$NumPortf = 30 # number of portfolios in efficient frontier
+Options$FrontierSpan = array( c(.1 , .8) ) # range of normalized exp.vals. spanned by efficient frontier
+Options$DeltaNeutral = 1 # 1 = portfolio is delta neutral, 0 = no constraint
+Options$Budget = 0 # initial capital
+Options$Limit = 10000 # limit to absolute value of each investment
+
+# The following commands loads the PnL object as calculated by the MATLAB option-pricing model
+ # This is for testing purposes so we can replicate results with the Meucci program
+load("C:\\Users\\Ram Ahluwalia\\Documents\\Applications\\R\\r-user\\MeucciButterflyPnL.rda")
+
+optimalPortfolios = LongShortMeanCVaRFrontier( PnL , p , Butterflies , Options )
+
+View( optimalPortfolios ) # Note that composition is measured in dollars. Here we are short GOOG_vol_91 and long GOOG_vol_182
+
+# plot efficient frontier ( Exp , CVaR , w )
+
+###################################################################
+# process views (this is the core of the Entropy Pooling approach
+###################################################################
+p_1 = ViewImpliedVol( X , p )
+ # View 1 (inequality view): bearish on on 2m-6pm implied volaility spread for Google
+ # Note that the mean-CVaR efficient frontier is long based on the reference model
+ # After processing the view, the G6m-G2m spread is now short in the mean-CVaR long-short efficient frontier
+
+p_2 = ViewRealizedVol( X , p )
+ # view 2 (relative inequality view on median): 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
+
+p_3 = ViewCurveSlope( X , p )
+ # view 3 (equality view - expectations and binding constraints): slope of the yield curve will increase by 5 bp
+
+# assign confidence to the views and pool opinions
+ # sum of weights should be equal 1
+ # .35 is the weight on the reference model
+ # .2 is the confidence on View 1; .25 is the confidence on View 2; .2 is the confidence on View 3
+ c = matrix( c( .35 , .2 , .25 , .2 ) , nrow = 1 )
+
+p_= cbind( p , p_1 , p_2 , p_3 ) %*% t(c) # compute the uncertainty weighted posterior probabilities
+
+###################################################################
+# compute mean-risk efficient frontier based on posterior market distribution
+###################################################################
+optimalPortfoliosPosterior = LongShortMeanCVaRFrontier( PnL , p_ , Butterflies , Options )
+
+View( optimalPortfoliosPosterior ) # Note that composition is measured in dollars. Now we are long GOOG_vol_91, and short Goog_vol_182
+
+# plot efficient frontier ( Exp_ , CVaR_ , w_ )
+PlotFrontier(Exp_,SDev_,w_)
+
+###################################################################
+# tests
+###################################################################
+p_3b = ViewCurveSlope( X , p )
+p_4 = ViewCurveSlopeTest( X , p )
+
+ViewCurveSlopeTest = function( X , p )
+ {
+ J = nrow( X ) ; K = ncol( X )
+
+ # constrain probabilities to sum to one...
+ Aeq = ones( 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 = .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_ )
+ }
Added: pkg/PerformanceAnalytics/sandbox/Meucci/MeucciCmaCopula.R
===================================================================
--- pkg/PerformanceAnalytics/sandbox/Meucci/MeucciCmaCopula.R (rev 0)
+++ pkg/PerformanceAnalytics/sandbox/Meucci/MeucciCmaCopula.R 2012-04-23 20:30:33 UTC (rev 1907)
@@ -0,0 +1,410 @@
+#' Risk and Asset Allocation functions
+#'
+#' \tabular{ll}{
+#' Package: \tab Atillio Meucci Risk functions\cr
+#' Type: \tab Package\cr
+#' Version: \tab 0.1\cr
+#' Date: \tab 2011-11-22\cr
+#' LazyLoad: \tab yes\cr
+#' }
+#'
+#' Contains a collection of functions to implement code from Atillio Meucci's text
+#' "Risk and Asset Allocation". See www.symmys.com for latest MATLAB source code and
+#' more information. Special thanks to Atillio Meucci and www.symmys.com
+#' for providing the MATLAB source on which most of the functions are based
+#'
+#' @name RiskAssetAllocation-package
+#' @aliases RiskAssetAllocation CMAcopula
+#' @docType package
+#' @title Implements the CMA copula based on Meucci (2011)
+#' @author Ram Ahluwalia \email{rahluwalia@@gmail.com}
+#' @references
+#' \url{http://www.symmys.com}
+#' @keywords package meucci copula entropy-pooling entropy CMA panic
+#' @examples
+#' print(" This is a test example at the top of the meuccicopula.R file ")
+# library(roxygen)
+library(roxygen2)
+roxygen:::roxygen()
+
+# Depends: \tab \cr matlab MASS roxygen2 trust pracma Hmisc Matrix \cr
+
+# fix documentation for @return in PanicCopula()
+# plot shape of the copula (perspective or surface plot)
+# fix MvnRnd function (Schur decomposition)
+# fix warnings
+
+#' Generates normal simulations whose sample moments match the population moments
+#'
+#' Adapted from file 'MvnRnd.m'. Most recent version of article and code available at http://www.symmys.com/node/162
+#' see A. Meucci - "Simulations with Exact Means and Covariances", Risk, July 2009
+#'
+#' @param M a numeric indicating the sample first moment of the distribution
+#' @param S a covariance matrix
+#' @param J a numeric indicating the number of trials
+#'
+#' @author Ram Ahluwalia \email{rahluwalia@@gmail.com}
+#' @references
+#' \url{http://www.symmys.com}
+#' TODO: Add Schur decomposition. Right now function is only sampling from mvrnorm so sample moments do no match population moments
+#' I have sample code commented out below to implement this correctly but I require a function that returns the unitaryMatrix from a Schur decomposition
+MvnRnd = function( M , S , J )
+ {
+ library(MASS)
+ X = MASS::mvrnorm( n = J , mu = M , Sigma = S ) # Todo: need to swap with Meucci function and Schur method
+ return( X = X )
+
+ # # compute sample covariance: NOTE defined as "cov(Y,1)", not as "cov(Y)"
+ # S_ = cov( Y , 1 )
+ #
+ # # solve Riccati equation using Schur method
+ # zerosMatrix = matrix( rep( 0 , length( N * N ) ) , nrow = N )
+ # # define the Hamiltonian matrix
+ # H1 = cbind( zerosMatrix , -1*S_ )
+ # H2 = cbind( -S , zerosMatrix )
+ # H = rbind( H1 , H2 )
+ # # perform its Schur decomposition.
+ # # TODO: check that the result returns real eigenvalues on the diagonal. ?Schur seems to give an example with real eigenvalues
+ # schurDecomp = Schur( H )
+ # T = SchurDecomp
+ # # U_ = unitaryMatrix??? TODO: Find a function in R that returns the unitaryMatrix from a Schur decomposition
+ # # U = ordschur(U_,T_,'lhp')
+ # # U_lu = U(1:N,1:N)
+ # # U_ld = U(N+1:end,1:N)
+ # # B = U_ld/U_lu
+ #
+ # # affine transformation to match mean and covariances
+ # # X = Y%*%B + repmat(M', J , 1 )
+ #
+ }
+
+
+#' CMA separation. Decomposes arbitrary joint distributions (scenario-probabilities) into their copula and marginals
+#'
+#' The CMA separation step attains from the cdf "F" for the marginal "X", the scenario-probabilities representation
+#' of the copula (cdf of U: "F") and the inter/extrapolation representation of the marginal CDF's
+#'
+#' Separation step of Copula-Marginal Algorithm (CMA)
+#' Meucci A., "New Breed of Copulas for Risk and Portfolio Management", Risk, September 2011
+#' Most recent version of article and code available at http://www.symmys.com/node/335
+#'
+#' @param X A matrix where each row corresponds to a scenario/sample from a joint distribution.
+#' Each column represents the value from a marginal distribution
+#' @param p A 1-column matrix of probabilities of the Jth-scenario joint distribution in X
+#'
+#' @return xdd a JxN matrix where each column consists of each marginal's generic x values in ascending order
+#' @return udd a JxN matrix containing the cumulative probability (cdf) for each marginal by column - it is rescaled by 'l' to be <1 at the far right of the distribution
+#' can interpret 'udd' as the probability weighted grade scenarios (see formula 11 in Meucci)
+#' @return U a copula (J x N matrix) - the joint distribution of grades defined by feeding the original variables X into their respective marginal CDF
+#'
+#'
+#' @author Ram Ahluwalia \email{rahluwalia@@gmail.com}
+#' @references
+#' \url{http://www.symmys.com}
+CMAseparation = function( X , p )
+{
+# @example test = cbind( seq( 102 , 1 , -2 ) , seq( 100 , 500 , 8 ) )
+# @example prob = rep(.02 , 51 )
+# @example CMAseparation( test , prob )
+
+library(pracma)
+# preprocess variables
+ X = as.matrix( X ) # J x N matrix. One column for each marginal distribution, J scenarios
+ p = as.matrix( p ) # J x 1 matrix of twisted probabilities
+ J = nrow( X )
+ N = ncol( X )
+ l = J / ( J + 1 ) # Laplace
+
+ # To ensure that no scenario in the prior is strictly impossible under the prior distribution we take the max of the prior or epsilon
+ # This is necessary otherwise one cannot apply Bayes rule if the prior probability of a scenario is 0
+ p = apply( cbind( p , 1 / J * 10^(-8) ) , 1 , max ) # return a J x 1 matrix where each row is the max of the row in p or 1 / J * 10^(-8)
+
+ p = p / sum(p) # Return a re-scaled vector of probabilities summing to 1
+ u = 0 * X # initialize a J x N matrix of 0's
+ U = 0 * X
+
+# begin core algorithm
+ # initialize empty J x N matrices
+ x = matrix( nrow = J , ncol = N )
+ Indx = matrix( nrow = J , ncol = N )
+
+ for ( n in 1:N ) # for each marginal...
+ {
+ x[ , n ] = sort( X[ , n ] , decreasing = FALSE ) # a JxN matrix where each column consists of each marginal's generic x values in ascending order
+ Indx[ , n ] = order( X[ , n ] ) # a JxN matrix. Each column contains the rank of the associated generic X in that column
+ }
+
+ for ( n in 1:N ) # for each marginal...
+ {
+ I = Indx[ , n ] # a Jx1 vector consisting of a given marginal's rank/index
+ cum_p = cumsum( p[I] ) # compute cdf. The head is a small number, and the tail value is 1
+ u[ , n ] = cum_p * l # 'u' will contain cdf for each marginal by column - it is rescaled by 'l' to be <1 at the far right of the distribution
+
+ # For each marginal distribution, CMA creates from each grid of scenario pairs { x-n,j , u-n,j } a function I{ x-n,j , u-n,j }
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/returnanalytics -r 1907
More information about the Returnanalytics-commits
mailing list