[Returnanalytics-commits] r2112 - pkg/PerformanceAnalytics/sandbox/Meucci/R

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

Author: mkshah
Date: 2012-07-05 06:20:17 +0200 (Thu, 05 Jul 2012)
New Revision: 2112

Modified:
   pkg/PerformanceAnalytics/sandbox/Meucci/R/HermiteGrid.R
Log:
Adding the kernel functions as defined by Meucci for the CaseStudy Demo

Modified: pkg/PerformanceAnalytics/sandbox/Meucci/R/HermiteGrid.R
===================================================================
--- pkg/PerformanceAnalytics/sandbox/Meucci/R/HermiteGrid.R	2012-07-04 23:12:29 UTC (rev 2111)
+++ pkg/PerformanceAnalytics/sandbox/Meucci/R/HermiteGrid.R	2012-07-05 04:20:17 UTC (rev 2112)
@@ -45,9 +45,7 @@
   if( n >= 2 ) {
     for (k in 2:n)
     {
-      firstPart = 2 *  cbind(t(p[k, 1:k]), 0)
-      secondPart = 2 * (k-1) * cbind(0, 0, t(p[k-1, 1:k-1]))
-      p[ k + 1, 1:(k + 1)] = firstPart  - secondPart
+      p[ k + 1, 1:(k + 1)] = 2 *  cbind(t(p[k, 1:k]), 0)  - 2 * (k-1) * cbind(0, 0, t(p[k-1, 1:k-1]))
     }
   }
     
@@ -61,4 +59,89 @@
   x = Re(x) 
                 
   return( x )
-}
\ No newline at end of file
+}
+
+# Silverman bandwidth
+kernelbw = function( xi )
+{
+  N      = length(xi)
+  prop   = 1.0
+  sig    = std( as.numeric( xi ) )
+  iqrSig = 0.7413 * IQR(xi, type = 5)
+
+  if (max(iqrSig) == 0) {
+    iqrSig = sig
+  }
+
+  bw = prop * min(sig, iqrSig) * N^(-1 / (4 + 1));
+  
+  return ( bw )
+}
+
+kernelcdf = function( x, xi, bw, wi )
+{
+  n = length(xi)
+  nargin <- length(as.list(match.call())) -1
+  if ( nargin < 4 || isempty( wi ) ) { wi = ones(n, 1) / n }
+
+  if ( nargin < 3 ) { bw = kernelbw(xi) }
+
+  p = zeros( size( x ) )
+  for ( i in 1:n ) {
+    p = p + exp( log( wi[i] ) + log( pnorm( x, xi[i], bw ) ) )
+  }
+  
+  return ( p )
+}
+
+kernelpdf = function( x, xi, bw, wi )
+{
+  n = length(xi)
+  nargin <- length(as.list(match.call())) -1
+  if ( nargin < 4 || isempty( wi ) ) { wi = ones(n, 1) / n }
+
+  if ( nargin < 3 ) { bw = kernelbw(xi) }
+
+  p = zeros( size( x ) )
+  for ( i in 1:n ) {
+    p = p + wi(i) * dnorm( x, xi(i), bw )
+  }
+}
+
+kernelinv = function( p, xi, bw, wi )
+{
+  nargin <- length(as.list(match.call())) -1
+  emptyMatrix = matrix( ,nrow = 0, ncol = 0)
+  if ( nargin < 4 || isempty(wi) ) { wi = emptyMatrix }
+  if ( nargin < 3 || isempty(bw) ) { bw = kernelbw(xi) }
+
+  sortp = sort(p)
+
+  if ( length(p) < 10 ) {
+    # case with only few points by treating each point seperately
+    x = zeros( dim( p ) )
+    for ( i in 1:length(p) ) {
+      x(i) = uniroot(function( x ) private_fun(x, xi, bw, wi, p[1]), 0)
+    }
+  }
+  else {
+  # case with many points by interpolation, find x_min and x_max
+  x_min = uniroot(function( x ) private_fun(x, xi, bw, wi, sortp[1]), 0 )
+  x_max = uniroot(function( x ) private_fun(x, xi, bw, wi, sortp[ nrow(sortp) ]), 0 )
+
+  # mesh for x values
+  x_ = seq( x_min - 0.1 * abs(x_min), x_max + 0.1 * abs(x_max), len = 500 )
+
+  # evaluates the mesh on these values
+  y_ = kernelcdf( x_, xi, bw, wi )
+
+  # interpolation
+  x = approx ( t( y_ ), t( x_ ), xout = t( p ) )
+  }
+}
+
+private_fun = function( x, xi, bw, wi, p )
+{
+  f = kernelcdf( x, xi, bw, wi ) - p
+  return( f )
+}            
\ No newline at end of file