Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
braverock
GitHub Repository: braverock/portfolioanalytics
 Path: blob/master/sandbox/FFEV/HermiteGrid_demo.R
1433 views
1
# This script compares the performance of plain Monte Carlo 

2
# versus grid in applying Entropy Pooling to process extreme views

3
# This script complements the article

4
#    "Fully Flexible Extreme Views" A. Meucci, D. Ardia, S. Keel available at www.ssrn.com

5
# The most recent version of this code is available at MATLAB Central - File Exchange

6
library(matlab)

7


8
####################################################################################

9
# Prior market model

10
####################################################################################

11
# analytical (normal) prior

12
emptyMatrix = matrix( nrow=0 , ncol=0 )

13
market.mu   = 0.0

14
market.sig2 = 1.0

15
market.pdf = function(x) dnorm( x , mean = market.mu , sd = sqrt(market.sig2) )

16
market.cdf = function(x) pnorm( x , mean = market.mu , sd = sqrt(market.sig2) )

17
market.rnd = function(x) rnorm( x , mean = market.mu , sd = sqrt(market.sig2) ) 

18
market.inv = function(x) qnorm( x , mean = market.mu , sd = sqrt(market.sig2) )

19


20
# numerical (Monte Carlo) prior 

21
monteCarlo = emptyMatrix

22
monteCarlo.J = 100000

23
monteCarlo.X = market.rnd( monteCarlo.J )

24
monteCarlo.p = normalizeProb( 1/monteCarlo.J * ones( monteCarlo.J , 1 ) )

25


26
# numerical (Gauss-Hermite grid) prior 

27
ghqMesh = emptyMatrix

28
load( "ghq1000" )

29


30
tmp = ( ghqx - min( ghqx ) ) / ( max( ghqx ) - min( ghqx ) ) # rescale GH zeros so they belong to [0,1]

31
epsilon = 1e-10

32
Lower = market.inv( epsilon )

33
Upper = market.inv( 1-epsilon )

34
ghqMesh.X  = Lower + tmp*(Upper-Lower) # rescale mesh

35


36
p = integrateSubIntervals(ghqMesh.X , market.cdf)

37
ghqMesh.p = normalizeProb(p)

38
ghqMesh.J = nrow(ghqMesh.X) 

39


40
####################################################################################

41
# Entropy posterior from extreme view on expectation

42
####################################################################################

43
# view of the analyst

44
view = emptyMatrix

45
view.mu = -3.0 

46


47
# analytical (known since normal model has analytical solution)

48
truePosterior = emptyMatrix 

49
truePosterior = Prior2Posterior( market.mu, 1, view.mu, market.sig2, 0 )

50
truePosterior$pdf = function(x) dnorm( x, truePosterior.mu , sqrt(truePosterior.sig2) )

51


52
# numerical (Monte Carlo)

53
Aeq = rbind( ones( 1 , monteCarlo.J ) , t(monteCarlo.X) )

54
beq = rbind( 1 , view.mu )

55
monteCarloOptimResult = EntropyProg( monteCarlo.p , emptyMatrix , emptyMatrix , Aeq , beq )

56


57
monteCarlo.p_ = monteCarloOptimResult$p_

58
monteCarlo.KLdiv = monteCarloOptimResult$optimizationPerformance$ml

59


60
# numerical (Gaussian-Hermite grid)

61
Aeq = rbind( ones( 1 , ghqMesh.J ) , t( ghqMesh.X ) )

62
beq = rbind( 1 , view.mu )

63
ghqMeshOptimResult = EntropyProg( ghqMesh.p , emptyMatrix , emptyMatrix , Aeq , beq )

64


65
ghqMesh.p_ = ghqMeshOptimResult$p_

66
ghqMesh.KLdiv = ghqMeshOptimResult$optimizationPerformance$ml

67


68
####################################################################################

69
# Plots

70
####################################################################################

71
xmin = min(ghqMesh.X)

72
xmax = max(ghqMesh.X)

73
ymax = 1.0

74
xmesh = t(linspace(xmin, xmax, ghqMesh.J))

75


76
# Monte Carlo

77
plotDataMC = pHist( monteCarlo.X , monteCarlo.p_ , 50 )

78
plot( plotDataMC$x , plotDataMC$f , type = "l" )

79


80
# Gauss Hermite Grid

81
plotDataGHQ = pHist(data.matrix(ghqMesh.X), ghqMesh.p_ , 50 )

82
plot( plotDataGHQ$x , plotDataGHQ$f , type = "l" )

83


84


85