[Returnanalytics-commits] r2374 - in pkg/Meucci: R demo
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Jun 19 13:31:59 CEST 2013
Author: xavierv
Date: 2013-06-19 13:31:59 +0200 (Wed, 19 Jun 2013)
New Revision: 2374
Added:
pkg/Meucci/R/LognormalCopulaPdf.R
pkg/Meucci/demo/S_AnalyzeLognormalCorrelation.R
pkg/Meucci/demo/S_AnalyzeNormalCorrelation.R
pkg/Meucci/demo/S_BivariateSample.R
pkg/Meucci/demo/S_DisplayLognormalCopulaPdf.R
Log:
- added S_DisplayLognormalCopulaPdf.R script and its required functions (fix)
Added: pkg/Meucci/R/LognormalCopulaPdf.R
===================================================================
--- pkg/Meucci/R/LognormalCopulaPdf.R (rev 0)
+++ pkg/Meucci/R/LognormalCopulaPdf.R 2013-06-19 11:31:59 UTC (rev 2374)
@@ -0,0 +1,35 @@
+
+#' Computes the pdf of the copula of the lognormal distribution at the generic point u in the unit hypercube,
+#' as described in A. Meucci, "Risk and Asset Allocation", Springer, 2005.
+#'
+#' @param u : [vector] (J x 1) grades
+#' @param Mu : [vector] (N x 1) location parameter
+#' @param Sigma : [matrix] (N x N) scatter parameter
+#'
+#' @return F_U : [vector] (J x 1) PDF values
+#'
+#' @references
+#' \url{http://}
+#' See Meucci's script for "LognormalCopulaPdf.m"
+#'
+#' @author Xavier Valls \email{flamejat@@gmail.com}
+#' @export
+
+LognormalCopulaPdf = function( u, Mu, Sigma )
+{
+ N = length( u );
+ s = sqrt( diag( Sigma ));
+
+ x = qlnorm( u, Mu, s );
+
+ Numerator = ( 2 * pi ) ^ ( -N / 2 ) * ( (det ( Sigma ) ) ^ ( -0.5 ) ) /
+ prod(x) * exp( -.5 * t(log(x) - Mu) %*% mldivide( Sigma , ( log( x ) - Mu ), pinv=FALSE ) );
+
+ fs = dlnorm( x, Mu, s);
+
+ Denominator = prod(fs);
+
+ F_U = Numerator / Denominator;
+
+ return ( F_U );
+}
\ No newline at end of file
Added: pkg/Meucci/demo/S_AnalyzeLognormalCorrelation.R
===================================================================
--- pkg/Meucci/demo/S_AnalyzeLognormalCorrelation.R (rev 0)
+++ pkg/Meucci/demo/S_AnalyzeLognormalCorrelation.R 2013-06-19 11:31:59 UTC (rev 2374)
@@ -0,0 +1,42 @@
+#' This script considers a bivariate lognormal market and display the correlation and the condition number of the
+#' covariance matrix, as described in A. Meucci, "Risk and Asset Allocation", Springer, 2005, Chapter 2.
+#'
+#' @references
+#' \url{http://}
+#' See Meucci's script for "S_AnalyzeLognormalCorrelation.m"
+#'
+#' @author Xavier Valls \email{flamejat@@gmail.com}
+#' @export
+
+
+###########################################################################################################################################
+### Set input parameters
+Mu = rbind( 0, 0 )
+s = c( 1, 1 );
+rhos = seq( -0.99, 0.99, 0.01 );
+nrhos = length( rhos );
+Cs = array( NaN, nrhos );
+CRs = array( NaN, nrhos );
+
+###########################################################################################################################################
+### Iterate of rho values and compute the correlation and condition number
+
+for ( n in 1 : nrhos )
+{
+ rho = rhos[ n ] ;
+ Sigma = rbind( c(s[1]^2, rho * s[1] * s[2]), c(rho * s[1] * s[2], s[2]^2) );
+
+ S = LognormalParam2Statistics(Mu, Sigma);
+
+ Lambda = eigen(S$Covariance);
+
+ Cs[ n ] = S$Correlation[1, 2];
+ CRs[ n ] = min(Lambda$values) / max(Lambda$values);
+}
+
+###########################################################################################################################################
+### Display the results
+
+par( mfrow = c( 2, 1) );
+plot( rhos, Cs, xlab = "r", ylab = "correlation", type ="l" );
+plot( rhos, CRs, xlab = "r", ylab = "condition ratio", type ="l" );
Added: pkg/Meucci/demo/S_AnalyzeNormalCorrelation.R
===================================================================
--- pkg/Meucci/demo/S_AnalyzeNormalCorrelation.R (rev 0)
+++ pkg/Meucci/demo/S_AnalyzeNormalCorrelation.R 2013-06-19 11:31:59 UTC (rev 2374)
@@ -0,0 +1,46 @@
+#' This script considers a bivariate normal market and display the correlation and the condition number of the
+#' covariance matrix, as described in A. Meucci, "Risk and Asset Allocation", Springer, 2005, Chapter 2.
+#'
+#' @references
+#' \url{http://}
+#' See Meucci's script for "S_AnalyzeNormalCorrelation.m"
+#'
+#' @author Xavier Valls \email{flamejat@@gmail.com}
+#' @export
+
+###################################################################################################################
+### Set input parameters
+
+Mu = rbind( 0, 0 )
+s = c( 1, 1 );
+
+rhos = seq( -0.99, 0.99, 0.01 );
+nrhos = length( rhos );
+
+Cs = array( NaN, nrhos );
+CRs = array( NaN, nrhos );
+
+
+###################################################################################################################
+### Iterate of rho values and compute the correlation and condition numberfor ( n in 1 : nrhos )
+
+for ( n in 1 : nrhos )
+{
+ rho = rhos[ n ] ;
+ Sigma = rbind( c(s[1]^2, rho * s[1] * s[2]), c(rho * s[1] * s[2], s[2]^2) );
+
+ Covariance = Sigma;
+ Standard_Deviation = sqrt( diag( Covariance ) );
+ Correlation = diag( 1 / Standard_Deviation ) %*% Covariance %*% diag( 1 / Standard_Deviation );
+
+ Lambda = eigen( Covariance );
+
+ Cs[n] = Correlation[ 1, 2 ];
+ CRs[n] = min( Lambda$values ) / max( Lambda$values );
+}
+
+###################################################################################################################
+### Display the results
+par( mfrow = c( 2, 1) );
+plot( rhos, Cs, xlab = "r", ylab = "correlation", type ="l" );
+plot( rhos, CRs, xlab = "r", ylab = "condition ratio", type ="l" );
Added: pkg/Meucci/demo/S_BivariateSample.R
===================================================================
--- pkg/Meucci/demo/S_BivariateSample.R (rev 0)
+++ pkg/Meucci/demo/S_BivariateSample.R 2013-06-19 11:31:59 UTC (rev 2374)
@@ -0,0 +1,117 @@
+#' This script generates draws from a bivariate distribution with different marginals,
+#' as described in A. Meucci, "Risk and Asset Allocation", Springer, 2005, Chapter 2.
+#'
+#' @references
+#' \url{http://}
+#' See Meucci's script for "S_AnalyzeLognormalCorrelation.m"
+#'
+#' @author Xavier Valls \email{flamejat@@gmail.com}
+#' @export
+
+library(mvtnorm);
+library(latticeExtra);
+
+###################################################################################################################
+### input parameters
+
+nSim = 10000;
+
+# input for bivariate normal distribution
+
+NormCorr = -0.8;
+NormStDev = rbind( 1, 3 ); # NOTE: this input plays no role in the final output
+NormExpVal = rbind( -2, 5 ); # NOTE: this input plays no role in the final output
+
+# input for first marginal
+nu_1 = 9;
+sigmasq_1 = 2;
+
+mu_2 = 0;
+sigmasq_2 = 0.04;
+
+# input for second marginal
+nu_2 = 7;
+
+###################################################################################################################
+### Generate draws from a bivariate normal distribution
+
+NormCorrMatrix = rbind( c( 1, NormCorr ), c( NormCorr, 1 ));
+NormCovMatrix = diag( c( NormStDev ) ) %*% NormCorrMatrix %*% diag( c( NormStDev) );
+
+Z = rvnorm( nSim, NormExpVal, NormCovMatrix );
+
+Z_1 = Z[, 1];
+Z_2 = Z[, 2];
+
+# display marginals: as expected, they are normal
+
+NumBins = round(10 * log(nSim));
+par( mfrow = c( 2, 1) );
+hist( Z_1, NumBins, xlab = "normal 1", ylab = "" );
+hist( Z_2, NumBins, xlab = "normal 2", ylab = "" );
+
+plot( Z_1, Z_2, type = "p", xlab = "normal 1", ylab = "normal 2" );
+
+
+# 3d histograms
+
+NumBins2D = round(sqrt(100 * log(nSim)));
+Z_3 = table( cut (Z_1, NumBins2D ), cut ( Z_2, cloud ));
+
+cloud( Z_3, panel.3d.cloud = panel.3dbars, scales = list( arrows = FALSE, just = "right" ),
+ xlab = "normal 1", ylab = "normal 2", zlab="", main = "pdf normal" );
+
+###################################################################################################################
+### Generate draws from the copula
+
+U_1 = pnorm( Z[ , 1 ], NormExpVal[ 1 ], NormStDev[ 1 ]); # grade 1
+U_2 = pnorm( Z[ , 2 ], NormExpVal[ 2 ], NormStDev[ 2 ]); # grade 2
+U = c( U_1, U_2 ); # joint realizations from the required copula
+
+# plot copula
+NumBins = round(10 * log(nSim));
+par( mfrow = c( 2, 1) );
+hist( U_1, NumBins, xlab = "grade 1", ylab = "", main = "" );
+hist( U_2, NumBins, xlab = "grade 2", ylab = "", main = "" );
+
+# joint sample
+plot(U_1, U_2, xlab="grade 1", ylab="grade 2" );
+
+# 3d histogram
+NumBins2D = round(sqrt(100 * log(nSim)));
+U_3 = table( cut (U_1, NumBins2D ), cut ( U_2, NumBins2D ));
+cloud( U_3, panel.3d.cloud = panel.3dbars, scales = list( arrows = FALSE, just = "right" ),
+ xlab = "grade 1", ylab = "grade 2", zlab="", main = "pdf copula" );
+
+###################################################################################################################
+### Generate draws from the joint distribution
+a = nu_1 / 2;
+b = 2 * sigmasq_1;
+X_1 = qgamma( U_1, a, b );
+
+sigma_2 = sqrt( sigmasq_2 );
+X_2 = qlnorm( U_2, mu_2, sigma_2 );
+
+X = C(X_1, X_2); # joint realizations from the required distribution
+
+###################################################################################################################
+### Plot joint distribution
+# marginals: as expected, the histograms (pdf's) do NOT change as NormCorr varies
+
+NumBins = round(10 * log(nSim));
+
+
+par( mfrow = c( 2, 1) );
+# Student t distribution
+hist( X_1, NumBins, xlab = "gamma", ylab = "", main = "" );
+# chi-square distribution
+hist( X_2, NumBins, xlab = "lognormal", ylab = "", main = "" );
+
+# joint sample
+plot(X_1, X_2, xlab="gamma", ylab="lognormal" );
+
+# 3d histogram
+NumBins2D = round(sqrt(100 * log(nSim)));
+X_3 = table( cut (X_1, NumBins2D ), cut ( X_2, NumBins2D ));
+cloud( X_3, panel.3d.cloud = panel.3dbars, scales = list( arrows = FALSE, just = "right" ),
+ xlab = "gamma", ylab = "lognormal", zlab="", main = "pdf joint distribution" );
\ No newline at end of file
Added: pkg/Meucci/demo/S_DisplayLognormalCopulaPdf.R
===================================================================
--- pkg/Meucci/demo/S_DisplayLognormalCopulaPdf.R (rev 0)
+++ pkg/Meucci/demo/S_DisplayLognormalCopulaPdf.R 2013-06-19 11:31:59 UTC (rev 2374)
@@ -0,0 +1,43 @@
+#'This script displays the pdf of the copula of a lognormal distribution, as described
+#' in A. Meucci, "Risk and Asset Allocation", Springer, 2005, Chapter 2.
+#'
+#' @references
+#' \url{http://}
+#' See Meucci's script for "S_DisplayLognormalCoulaPdf.m"
+#'
+#' @author Xavier Valls \email{flamejat@@gmail.com}
+#' @export
+
+#############################################################################################################
+### Input parameters
+Mu = rbind( 100.0, -30.0 );
+r = 0.8;
+sigmas = rbind( 1000, 0.01 );
+nu = 100;
+Sigma = diag( c( sigmas ) ) %*% rbind( c( 1, r ), c( r, 1 ) ) %*% diag( c( sigmas ) );
+
+#############################################################################################################
+### Grid
+GridSide1 = seq( 0.05, 0.95, 0.05 );
+GridSide2 = GridSide1;
+nMesh = length(GridSide1);
+
+#############################################################################################################
+### Compute pdf of copula
+
+f_U = matrix( NaN, nMesh, nMesh);
+
+for ( j in 1 : nMesh )
+{
+ for ( k in 1 : nMesh)
+ {
+ u = c( GridSide1[ j ], GridSide2[ k ] );
+ f_U[ j, k ] = LognormalCopulaPdf(u, Mu, Sigma);
+ }
+}
+
+#mesh representation
+
+persp(GridSide1,GridSide2, f_U,
+ theta = 7 * 45, phi = 30, expand=0.6, col='lightblue', shade=0.75, ltheta=120,
+ ticktype='detailed', xlab = "U_1", ylab = "U_2", zlab = "copula pdf" );
More information about the Returnanalytics-commits
mailing list