[Returnanalytics-commits] r2494 - in pkg/Meucci: . R demo man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Jul 3 08:38:24 CEST 2013
Author: xavierv
Date: 2013-07-03 08:38:24 +0200 (Wed, 03 Jul 2013)
New Revision: 2494
Added:
pkg/Meucci/R/BlackScholesCallPrice.R
pkg/Meucci/R/InterExtrapolate.R
pkg/Meucci/demo/S_CallsProjectionPricing.R
pkg/Meucci/man/BlackScholesCallPrice.Rd
pkg/Meucci/man/InterExtrapolate.Rd
Modified:
pkg/Meucci/DESCRIPTION
pkg/Meucci/NAMESPACE
Log:
-added S_CallsProjectionPricing demo script and its associated functions
Modified: pkg/Meucci/DESCRIPTION
===================================================================
--- pkg/Meucci/DESCRIPTION 2013-07-03 01:55:18 UTC (rev 2493)
+++ pkg/Meucci/DESCRIPTION 2013-07-03 06:38:24 UTC (rev 2494)
@@ -73,3 +73,5 @@
'TwoDimEllipsoid.R'
'PerformIidAnalysis.R'
'SimulateJumpDiffusionMerton.R'
+ 'BlackScholesCallPrice.R'
+ 'InterExtrapolate.R'
Modified: pkg/Meucci/NAMESPACE
===================================================================
--- pkg/Meucci/NAMESPACE 2013-07-03 01:55:18 UTC (rev 2493)
+++ pkg/Meucci/NAMESPACE 2013-07-03 06:38:24 UTC (rev 2494)
@@ -1,3 +1,4 @@
+export(BlackScholesCallPrice)
export(Central2Raw)
export(CMAcombination)
export(CMAseparation)
@@ -11,6 +12,7 @@
export(GenerateLogNormalDistribution)
export(hermitePolynomial)
export(integrateSubIntervals)
+export(InterExtrapolate)
export(linreturn)
export(LognormalCopulaPdf)
export(LognormalMoments2Parameters)
Added: pkg/Meucci/R/BlackScholesCallPrice.R
===================================================================
--- pkg/Meucci/R/BlackScholesCallPrice.R (rev 0)
+++ pkg/Meucci/R/BlackScholesCallPrice.R 2013-07-03 06:38:24 UTC (rev 2494)
@@ -0,0 +1,33 @@
+#' Compute the Black-Scholes price of a European call option
+#' as described in A. Meucci, "Risk and Asset Allocation", Springer, 2005.
+#'
+#' @param spot : [scalar] spot price of underlying
+#' @param K : [scalar] strike of the call optioon
+#' @param r : [scalar] risk free rate as a fraction
+#' @param vol : [scalar] volatility of the underlying as a fraction
+#' @param T : [scalar] time to maturity in years
+#'
+#' @return c : [scalar] price of European call(s)
+#' @return delta : [scalar] delta of the call(s)
+#' @return cash : [scalar] cash held in a replicating portfolio
+#'
+#' @note
+#' Code is vectorized, so the inputs can be vectors or matrices (but sizes must match)
+#'
+#' @references
+#' \url{http://symmys.com/node/170}
+#' See Meucci's script for "BlackScholesCallPrice.m"
+#'
+#' @author Xavier Valls \email{flamejat@@gmail.com}
+#' @export
+
+BlackScholesCallPrice = function(spot, K, r, vol, T)
+{
+ d1 = ( log( spot / K ) + ( r + vol * vol / 2) * T) / (vol * sqrt(T));
+ d2 = d1 - vol * sqrt(T);
+ delta = pnorm(d1);
+ cash = -K * exp( -r * T ) * pnorm( d2 );
+ c = spot * delta + cash;
+
+ return( list( c = c, delta = delta, cash = cash ) );
+}
\ No newline at end of file
Added: pkg/Meucci/R/InterExtrapolate.R
===================================================================
--- pkg/Meucci/R/InterExtrapolate.R (rev 0)
+++ pkg/Meucci/R/InterExtrapolate.R 2013-07-03 06:38:24 UTC (rev 2494)
@@ -0,0 +1,174 @@
+#' Interpolate and extrapolate using n-linear interpolation (tensor product linear).
+#'
+#' @param V : [array] p-dimensional array to be interpolated/extrapolated at the list of points in the array Xi.
+# interpne will work in any number of dimensions >= 1
+#' @param Xi : [array] (n x p) array of n points to interpolate/extrapolate. Each point is one row of the array Xi.
+#' @param nodelist : [cell array] (optional) cell array of nodes in each dimension.
+# If nodelist is not provided, then by default I will assume nodelist[[i]] = 1:size(V,i). The nodes in
+# nodelist need not be uniformly spaced.
+#' @param method : [string] (optional) chacter string, denotes the interpolation method used. default method = 'linear'
+# 'linear' --> n-d linear tensor product interpolation/extrapolation
+# 'nearest' --> n-d nearest neighbor interpolation/extrapolation
+# in 2-d, 'linear' is equivalent to a bilinear interpolant
+# in 3-d, it is commonly known as trilinear interpolation.
+#'
+#' @return Vpred : [array] (n x 1) array of interpolated/extrapolated values
+#'
+#' @note
+#' Initially written by John D'Errico
+#' Vpred = interpne(V,Xi)
+#' Vpred = interpne(V,Xi,nodelist)
+#' Vpred = interpne(V,Xi,nodelist,method)
+#' Extrapolating long distances outside the support of V is rarely advisable.
+#'
+#' @examples
+#'
+#' [x1,x2] = meshgrid(0:.2:1);
+#' z = exp(x1+x2);
+#' Xi = rand(100,2)*2-.5;
+#' Zi = interpne(z,Xi,{0:.2:1, 0:.2:1},'linear');
+#' surf(0:.2:1,0:.2:1,z)
+#' hold on
+#' plot3(Xi(:,1),Xi(:,2),Zi,'ro')
+#'
+#' @references
+#' \url{http://symmys.com/node/170}
+#' See Meucci's script for "InterExtrapolate.R"
+#'
+#' @author Xavier Valls \email{flamejat@@gmail.com}
+#' @export
+
+InterExtrapolate = function( V, Xi, nodelist, method, ...)
+{
+ # get some sizes
+
+ vsize = dim( V );
+ ndims = length( vsize );
+ nargin = length( match.call() ) -1;
+ if( ndims != ncol( Xi ) ) stop("Xi is not compatible in size with the array V for interpolation.")
+
+ # default for nodelist
+ if ( nargin < 3 || length(nodelist) == 0 )
+ {
+ nodelist= vector( "list", ndims );
+
+ for( i in 1 : ndims )
+ {
+ nodelist[[ i ]] = rbind( 1 : vsize[ i ] );
+ }
+ }
+
+ if ( length( nodelist ) != ndims )
+ stop( "nodelist is incompatible with the size of V.")
+
+ nll = lapply( nodelist, length );
+
+ if( any( nll!= vsize) )
+ stop( "nodelist is incompatible with the size of V." )
+
+ # get deltax for the node spacing
+ dx = nodelist;
+
+ for( i in 1 : ndims )
+ {
+ nodelist[[i]] = nodelist[[i]][]; # not sure about this doing anything
+ dx[[i]] = diff(nodelist[[i]][,]);
+ if ( any( dx[[i]] <= 0) )
+ stop( "The nodes in nodelist must be monotone increasing." );
+ }
+
+ # check for method
+ if ( nargin < 4 ) method = 'linear';
+
+
+
+ if( ! is.character(method) ) stop("method must be a character string if supplied.");
+
+
+ validmethod = c( "linear", "nearest");
+
+ if(!any(validmethod == method ) )
+ stop(paste(" No match found for method = ", method))
+
+ # Which cell of the array does each point lie in?
+ # This includes extrapolated points, which are also taken
+ # to fall in a cell. histc will do all the real work.
+
+ ind = matrix( 0, nrow(Xi), ndims);
+
+ for( i in 1 : ndims)
+ {
+ hc = histc(Xi[ , i ], nodelist[[i]]); ##ok<ASGLU>
+
+ # catch any point along the very top edge.
+ hc$bin[ hc$bin == vsize[ i ] ] = vsize[ i ] - 1;
+
+ ind[ , i ] = hc$bin;
+
+ k = which( hc$bin == 0);
+
+ # look for any points external to the nodes
+ if( !(length(k)==0) )
+ {
+ # bottom end
+ ind[ k[ Xi[ k, i] < nodelist[[i]][ 1 ]], i ] = 1;
+
+ # top end
+ ind[ k[ Xi[ k, i] > nodelist[[i]][ length(nodelist[[1]][1]) ]], i ] = vsize[ i ] - 1;
+ }
+ }
+
+ # where in each cell does each point fall?
+ t = matrix( 0, nrow(Xi), ndims);
+
+ for( i in 1 : ndims)
+ {
+ t[ , i ] = (Xi[ , i ] - nodelist[[i]][ ind[ , i ] ] )/ dx[[i]][ind[ , i ]];
+ }
+
+ sub = cumprod( c( 1 ,vsize[ 1 : ( length( vsize ) - 1 ) ] ) );
+ base = 1 + ( ind-1 ) %*% sub;
+
+ # which interpolation method do we use?
+
+ switch( method,
+ nearest =
+ {
+ # nearest neighbor is really simple to do.
+ t = round(t);
+ t[ t > 1 ] = 1;
+ t[ t < 0 ] = 0;
+
+ Vpred = V[ base + t %*% sub ];
+ },
+
+ linear =
+ {
+ # tensor product linear is not too nasty.
+ Vpred = matrix( 0, nrow(Xi), 1);
+ # define the 2^ndims corners of a hypercube (MATLAB's corners = (dec2bin(0:(2^ndims-1))== '1');)
+ corners = lapply( strsplit( intToBin ( 0 : ( 2^ndims - 1 ) ), split=""), as.integer );
+
+ nc = length( corners );
+
+ for( i in 1 : nc )
+ {
+ #accessing
+ s = V[ base + (corners[[i]] %*% sub)[1]];
+ for( j in 1 : ndims )
+ {
+ # this will work for extrapolation too
+ if( corners[[i]][ j ] == 0 ){
+ s = s * ( 1 - t[ , j ] );
+ }else
+ {
+ s = s * t[ , j ];
+ }
+ }
+
+ Vpred = Vpred + s;
+ }
+ } ) # end switch method
+
+ return( Vpred );
+}
\ No newline at end of file
Added: pkg/Meucci/demo/S_CallsProjectionPricing.R
===================================================================
--- pkg/Meucci/demo/S_CallsProjectionPricing.R (rev 0)
+++ pkg/Meucci/demo/S_CallsProjectionPricing.R 2013-07-03 06:38:24 UTC (rev 2494)
@@ -0,0 +1,105 @@
+library(mvtnorm);
+library(pracma);
+
+#'This script projects the distribution of the market invariants for the derivatives market
+#'Then it computes the distribution of prices at the investment horizon as described in A. Meucci,
+#'"Risk and Asset Allocation", Springer, 2005, Chapter 3.
+#'
+#' @references
+#' \url{http://symmys.com/node/170}
+#' See Meucci's script for "S_CallsProjectionPricing.m"
+#'
+#' @author Xavier Valls \email{flamejat@@gmail.com}
+
+
+##################################################################################################################
+### Load data
+
+# load 'spot' for underlying and current vol surface, given by
+# 'impVol' for different 'days2Maturity' and 'moneyness' (K/S)
+load("../data/implVol.Rda");
+
+##################################################################################################################
+### Inputs
+
+tau_tilde = 5; # estimation step (days)
+tau = 40; # time to horizon (days)
+Time2Mats = c( 100, 150, 200, 250, 300 ); # current time to maturity of call options in days
+Strikes = c( 850, 880, 910, 940, 970 ); # strikes of call options, same dimension as Time2Mat
+
+r_free = 0.04; # risk-free rate
+J = 10000; # number of simulations
+
+##################################################################################################################
+numCalls = length( Time2Mats );
+timeLength = length( implVol$spot );
+numSurfPoints = length( implVol$days2Maturity ) * length( implVol$moneyness );
+
+##################################################################################################################
+### Estimate invariant distribution assuming normality
+# variables in X are changes in log(spot) and changes in log(imp.vol)
+# evaluated at the 'numSurfPoints' points on the vol surface (vectorized).
+X = matrix( 0, timeLength - 1, numSurfPoints + 1 );
+# log-changes of underlying spot
+X[ , 1 ] = diff( log( implVol$spot ) );
+
+# log-changes of implied vol for different maturities
+impVolSeries = matrix( implVol$impVol, timeLength, numSurfPoints );
+for( i in 1 : numSurfPoints )
+{
+ X[ , i+1 ] = diff( log( impVolSeries[ , i ] ) );
+}
+
+muX = apply( X , 2, mean );
+SigmaX = cov( X );
+
+##################################################################################################################
+### Project distribution to investment horizon
+muX = muX * tau / tau_tilde;
+SigmaX = SigmaX * tau / tau_tilde;
+
+##################################################################################################################
+### Linearly interpolate the vol surface at the current time to obtain implied vol for the given calls today, and price the calls
+spot_T = implVol$spot[ length(implVol$spot) ];
+volSurf_T = drop( implVol$impVol[ length(implVol$impVol[, 1, 1] ), , ]);
+time2Mat_T = Time2Mats;
+moneyness_T = Strikes/spot_T;
+
+impVol_T = t( InterExtrapolate( volSurf_T,t( rbind( time2Mat_T, moneyness_T )), list( implVol$days2Maturity,implVol$moneyness ) ) );
+callPrice_T = BlackScholesCallPrice( spot_T, Strikes, r_free, impVol_T, Time2Mats/252 )$c;
+
+##################################################################################################################
+### Generate simulations at horizon
+X_ = rmvnorm( J, muX, SigmaX );
+
+# interpolate vol surface at horizon for the given calls
+spot_ = spot_T * exp(X_[ , 1 ] );
+impVol_ = matrix( 0, J, numCalls);
+for( j in 1 : J )
+{
+ volSurf = volSurf_T * exp( matrix(X_[ j, -1 ],length( implVol$days2Maturity ),length( implVol$moneyness )));
+ time2Mat_ = Time2Mats - tau;
+ moneyness_ = Strikes / spot_[ j ];
+ impVol_[ j, ] = t( InterExtrapolate( volSurf,t( rbind( time2Mat_T, moneyness_T )), list( implVol$days2Maturity,implVol$moneyness ) ));
+}
+
+# price the calls at the horizon
+callPrice_ = matrix( 0, J, numCalls );
+for( i in 1 : numCalls )
+{
+ callPrice_[ , i ] = BlackScholesCallPrice( spot_, Strikes[ i ], r_free, impVol_[ , i ], time2Mat_[ i ] / 252 )$c;
+}
+
+m = nrow( callPrice_ );
+n = ncol( callPrice_ );
+LinearRets = callPrice_ /kronecker( matrix( 1, J, 1), callPrice_T)-1
+
+NumBins = round(10 * log(J));
+
+for( i in 1 : numCalls)
+{
+ dev.new();
+ par( mfrow = c( 2 , 1));
+ hist( callPrice_[ , i ], NumBins, xlab = "call price");
+ plot(spot_, callPrice_[ ,i ], xlab = "spot price", ylab = "call price" );
+}
\ No newline at end of file
Added: pkg/Meucci/man/BlackScholesCallPrice.Rd
===================================================================
--- pkg/Meucci/man/BlackScholesCallPrice.Rd (rev 0)
+++ pkg/Meucci/man/BlackScholesCallPrice.Rd 2013-07-03 06:38:24 UTC (rev 2494)
@@ -0,0 +1,43 @@
+\name{BlackScholesCallPrice}
+\alias{BlackScholesCallPrice}
+\title{Compute the Black-Scholes price of a European call option
+as described in A. Meucci, "Risk and Asset Allocation", Springer, 2005.}
+\usage{
+ BlackScholesCallPrice(spot, K, r, vol, T)
+}
+\arguments{
+ \item{spot}{: [scalar] spot price of underlying}
+
+ \item{K}{: [scalar] strike of the call optioon}
+
+ \item{r}{: [scalar] risk free rate as a fraction}
+
+ \item{vol}{: [scalar] volatility of the underlying as a
+ fraction}
+
+ \item{T}{: [scalar] time to maturity in years}
+}
+\value{
+ c : [scalar] price of European call(s)
+
+ delta : [scalar] delta of the call(s)
+
+ cash : [scalar] cash held in a replicating portfolio
+}
+\description{
+ Compute the Black-Scholes price of a European call option
+ as described in A. Meucci, "Risk and Asset Allocation",
+ Springer, 2005.
+}
+\note{
+ Code is vectorized, so the inputs can be vectors or
+ matrices (but sizes must match)
+}
+\author{
+ Xavier Valls \email{flamejat at gmail.com}
+}
+\references{
+ \url{http://symmys.com/node/170} See Meucci's script for
+ "BlackScholesCallPrice.m"
+}
+
Added: pkg/Meucci/man/InterExtrapolate.Rd
===================================================================
--- pkg/Meucci/man/InterExtrapolate.Rd (rev 0)
+++ pkg/Meucci/man/InterExtrapolate.Rd 2013-07-03 06:38:24 UTC (rev 2494)
@@ -0,0 +1,53 @@
+\name{InterExtrapolate}
+\alias{InterExtrapolate}
+\title{Interpolate and extrapolate using n-linear interpolation (tensor product linear).}
+\usage{
+ InterExtrapolate(V, Xi, nodelist, method, ...)
+}
+\arguments{
+ \item{V}{: [array] p-dimensional array to be
+ interpolated/extrapolated at the list of points in the
+ array Xi.}
+
+ \item{Xi}{: [array] (n x p) array of n points to
+ interpolate/extrapolate. Each point is one row of the
+ array Xi.}
+
+ \item{nodelist}{: [cell array] (optional) cell array of
+ nodes in each dimension.}
+
+ \item{method}{: [string] (optional) chacter string,
+ denotes the interpolation method used. default method =
+ 'linear'}
+}
+\value{
+ Vpred : [array] (n x 1) array of
+ interpolated/extrapolated values
+}
+\description{
+ Interpolate and extrapolate using n-linear interpolation
+ (tensor product linear).
+}
+\note{
+ Initially written by John D'Errico Vpred = interpne(V,Xi)
+ Vpred = interpne(V,Xi,nodelist) Vpred =
+ interpne(V,Xi,nodelist,method) Extrapolating long
+ distances outside the support of V is rarely advisable.
+}
+\examples{
+[x1,x2] = meshgrid(0:.2:1);
+ z = exp(x1+x2);
+ Xi = rand(100,2)*2-.5;
+ Zi = interpne(z,Xi,{0:.2:1, 0:.2:1},'linear');
+ surf(0:.2:1,0:.2:1,z)
+ hold on
+ plot3(Xi(:,1),Xi(:,2),Zi,'ro')
+}
+\author{
+ Xavier Valls \email{flamejat at gmail.com}
+}
+\references{
+ \url{http://symmys.com/node/170} See Meucci's script for
+ "InterExtrapolate.R"
+}
+
More information about the Returnanalytics-commits
mailing list