[Returnanalytics-commits] r2121 - pkg/PerformanceAnalytics/sandbox/Meucci/demo

[noreply at r-forge.r-project.org]

Author: mkshah
Date: 2012-07-07 20:25:43 +0200 (Sat, 07 Jul 2012)
New Revision: 2121

Modified:
   pkg/PerformanceAnalytics/sandbox/Meucci/demo/HermiteGrid_CaseStudy.R
Log:
Completed and checked the results of the demo with Meucci's Matlab code

Modified: pkg/PerformanceAnalytics/sandbox/Meucci/demo/HermiteGrid_CaseStudy.R
===================================================================
--- pkg/PerformanceAnalytics/sandbox/Meucci/demo/HermiteGrid_CaseStudy.R	2012-07-07 05:34:55 UTC (rev 2120)
+++ pkg/PerformanceAnalytics/sandbox/Meucci/demo/HermiteGrid_CaseStudy.R	2012-07-07 18:25:43 UTC (rev 2121)
@@ -15,9 +15,9 @@
 # bandwidth
 bw = kernelbw(xi)
 
-% weights
+# weights
 lambda = log(2) / (n / 2)
-wi = as.matrix( exp( -lambda * ( n - seq( from=n, to=1 ) ) )
+wi = as.matrix( exp( -lambda * ( n - seq( from=n, to=1 ) ) ) )
 wi = as.matrix( apply( wi, 2, rev ) / sum( wi ) )
                 
 # Prior market model
@@ -27,8 +27,8 @@
 market.pdf = function ( x ) kernelpdf( x, xi, bw, wi )
 market.cdf = function ( x ) kernelcdf( x, xi, bw, wi )
 market.inv = function ( x ) kernelinv( x, xi, bw, wi )
-market.VaR95 = market.inv(0.05)
-market.CVaR95 = integrate( function( x ) (x %*% market.pdf( x ) ), -100, market.VaR95 ) / 0.05;
+market.VaR95 = market.inv( c(0.05) )
+market.CVaR95 = integrate( function( x ) x * market.pdf( x ), -100, market.VaR95 )$val / 0.05
 
 # numerical (Gauss-Hermite grid) prior 
 tmp = ( ghqx - min( ghqx ) )/( max( ghqx ) - min( ghqx ) ) # rescale GH zeros so they belong to [0,1]
@@ -39,9 +39,57 @@
                 
 p = integrateSubIntervals( X , market.cdf )
 p = normalizeProb( p )
-J = length( X )
+J = nrow( X )
                 
 # Entropy posterior from extreme view on mean and CVaR
+view.mu   = mean( xi ) - 1.0
+view.CVaR95 = market.CVaR95 - 1.0
                 
-                
-                
\ No newline at end of file
+# Netwton Raphson
+emptyMatrix = matrix( ,nrow = 0, ncol = 0)
+s = emptyMatrix
+idx = as.matrix( cumsum(p) <= 0.05 )
+s[1] = sum(idx)
+posteriorEntropy = EntropyProg(p, t( idx ), as.matrix( 0.05 ), rbind( rep(1, J), t( X ), t( idx * X ) ), rbind( 1, view.mu, 0.05 * view.CVaR95) )
+KLdiv = as.matrix( posteriorEntropy$optimizationPerformance$ml )
+p_ = posteriorEntropy$p_
+
+doStop = 0
+i = 1
+while ( !doStop ) {
+  i = i + 1
+
+  idx = cbind( matrix(1, 1, s[i - 1] ), matrix(0, 1, J - s[i-1] ) )
+  posteriorEntropy1 = EntropyProg(p, idx, as.matrix( 0.05 ), rbind( matrix(1, 1, J), t( X ), t( t(idx) * X) ), rbind( 1, view.mu, 0.05 * view.CVaR95 ) )
+  # [dummy, KLdiv_s]  = optimizeEntropy(p, [idx'; (idx .* X)'], [0.05; 0.05 * view.CVaR95], [ones(1, J); X'], [1; view.mu]);
+   
+  idx = cbind( matrix(1, 1, s[i - 1] + 1 ), matrix(0, 1, J - s[i - 1] - 1) )
+  posteriorEntropy2 = EntropyProg(p, idx, as.matrix( 0.05 ), rbind( matrix(1, 1, J), t( X ), t( t(idx) * X) ), rbind( 1, view.mu, 0.05 * view.CVaR95 ) )
+  # [dummy, KLdiv_s1] = optimizeEntropy(p, [idx'; (idx .* X)'], [0.05; 0.05 * view.CVaR95], [ones(1, J); X'], [1; view.mu]);
+   
+  idx = cbind( matrix(1, 1, s[i - 1] + 2 ), matrix(0, 1, J - s[i - 1] - 2) )
+  posteriorEntropy3 = EntropyProg(p, idx, as.matrix( 0.05 ), rbind( matrix(1, 1, J), t( X ), t( t(idx) * X) ), rbind( 1, view.mu, 0.05 * view.CVaR95 ) )
+  # [dummy, KLdiv_s2] = optimizeEntropy(p, [idx'; (idx .* X)'], [0.05; 0.05 * view.CVaR95], [ones(1, J); X'], [1; view.mu]);
+
+  # first difference
+  DE  = posteriorEntropy2$optimizationPerformance$ml - posteriorEntropy1$optimizationPerformance$ml
+  
+  # second difference
+  D2E = posteriorEntropy3$optimizationPerformance$ml - 2 * posteriorEntropy2$optimizationPerformance$ml + posteriorEntropy1$optimizationPerformance$ml
+  
+  # optimal s
+  s = rbind( s, round( s[i - 1] - (DE / D2E) ) ) 
+                                                                                             
+  tmp = emptyMatrix
+  idx  = cbind( matrix( 1, 1, s[i] ), matrix( 0, 1, J - s[i] ) )
+  tempEntropy = EntropyProg(p, idx, as.matrix( 0.05 ), rbind( matrix(1, 1, J), t( X ), t( t(idx) * X) ), rbind( 1, view.mu, 0.05 * view.CVaR95 ) )
+  # [tmp.p_, tmp.KLdiv] = optimizeEntropy(p, [idx'; (idx .* X)'], [0.05; 0.05 * view.CVaR95], [ones(1, J); X'], [1; view.mu]);
+  p_ = cbind( p_, tempEntropy$p_ )
+  KLdiv = rbind( KLdiv, tempEntropy$optimizationPerformance$ml ) #ok<*AGROW>
+
+  # if change in KLdiv less than one percent, stop
+  if( abs( ( KLdiv[i] - KLdiv[i - 1] ) / KLdiv[i - 1] )  < 0.01 ) { doStop = 1 }
+}                
+
+plot( t(X), p )
+plot( t(X), p_[,ncol(p_)])
\ No newline at end of file