-
Notifications
You must be signed in to change notification settings - Fork 1
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Merge pull request #3 from louisaslett/tiny-fixes
Small fixes in documentation and pkgdown, plus some tiny bug fixes and a big bug fix in opre_l()
- Loading branch information
Showing
15 changed files
with
252 additions
and
169 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -3,7 +3,7 @@ | |
|
|
||
| #' Financial options using a Milstein discretisation | ||
| #' | ||
| #' Financial options based on scalar geometric Brownian motion, similar to Mike Giles' MCQMC06 paper, using a Milstein discretisation. | ||
| #' Financial options based on scalar geometric Brownian motion, similar to Mike Giles' MCQMC06 paper, Giles (2008), using a Milstein discretisation. | ||
| #' | ||
| #' This function is based on GPL-2 C++ code by Mike Giles. | ||
| #' | ||
|
|
@@ -26,7 +26,7 @@ | |
| #' @author Mike Giles <[email protected]> | ||
| #' | ||
| #' @references | ||
| #' M.B. Giles. 'Improved multilevel Monte Carlo convergence using the Milstein scheme', p.343-358 in \emph{Monte Carlo and Quasi-Monte Carlo Methods 2006}, Springer, 2007. | ||
| #' Giles, M. (2008) 'Improved Multilevel Monte Carlo Convergence using the Milstein Scheme', in A. Keller, S. Heinrich, and H. Niederreiter (eds) \emph{Monte Carlo and Quasi-Monte Carlo Methods 2006}. Berlin, Heidelberg: Springer, pp. 343–358. Available at: \url{https://doi.org/10.1007/978-3-540-74496-2_20}. | ||
| #' | ||
| #' @examples | ||
| #' \dontrun{ | ||
|
|
@@ -39,49 +39,48 @@ | |
| #' # -- different random number generation | ||
| #' # -- switch to S_0=100 | ||
| #' | ||
| #' M <- 2 # refinement cost factor | ||
| #' N0 <- 200 # initial samples on coarse levels | ||
| #' Lmin <- 2 # minimum refinement level | ||
| #' Lmax <- 10 # maximum refinement level | ||
| #' | ||
| #' test.res <- list() | ||
| #' for(option in 1:5) { | ||
| #' if(option==1) { | ||
| #' if(option == 1) { | ||
| #' cat("\n ---- Computing European call ---- \n") | ||
| #' N <- 20000 # samples for convergence tests | ||
| #' L <- 8 # levels for convergence tests | ||
| #' Eps <- c(0.005, 0.01, 0.02, 0.05, 0.1) | ||
| #' } else if(option==2) { | ||
| #' } else if(option == 2) { | ||
| #' cat("\n ---- Computing Asian call ---- \n") | ||
| #' N <- 20000 # samples for convergence tests | ||
| #' L <- 8 # levels for convergence tests | ||
| #' Eps <- c(0.005, 0.01, 0.02, 0.05, 0.1) | ||
| #' } else if(option==3) { | ||
| #' } else if(option == 3) { | ||
| #' cat("\n ---- Computing lookback call ---- \n") | ||
| #' N <- 20000 # samples for convergence tests | ||
| #' L <- 10 # levels for convergence tests | ||
| #' Eps <- c(0.005, 0.01, 0.02, 0.05, 0.1) | ||
| #' } else if(option==4) { | ||
| #' } else if(option == 4) { | ||
| #' cat("\n ---- Computing digital call ---- \n") | ||
| #' N <- 200000 # samples for convergence tests | ||
| #' L <- 8 # levels for convergence tests | ||
| #' Eps <- c(0.01, 0.02, 0.05, 0.1, 0.2) | ||
| #' } else if(option==5) { | ||
| #' } else if(option == 5) { | ||
| #' cat("\n ---- Computing barrier call ---- \n") | ||
| #' N <- 200000 # samples for convergence tests | ||
| #' L <- 8 # levels for convergence tests | ||
| #' Eps <- c(0.005, 0.01, 0.02, 0.05, 0.1) | ||
| #' } | ||
| #' | ||
| #' test.res[[option]] <- mlmc.test(mcqmc06_l, M, N, L, N0, Eps, Lmin, Lmax, option=option) | ||
| #' test.res[[option]] <- mlmc.test(mcqmc06_l, N, L, N0, Eps, Lmin, Lmax, option = option) | ||
| #' | ||
| #' # plot results | ||
| #' plot(test.res[[option]]) | ||
| #' } | ||
| #' } | ||
| #' | ||
| #' # The level sampler can be called directly to retrieve the relevant level sums: | ||
| #' mcqmc06_l(l=7, N=10, option=1) | ||
| #' mcqmc06_l(l = 7, N = 10, option = 1) | ||
| #' | ||
| #' @export | ||
| mcqmc06_l <- function(l, N, option) { | ||
|
|
||
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -10,8 +10,9 @@ | |
| #' Consider a sequence \eqn{P_0, P_1, \ldots}, which approximates \eqn{P_L} with increasing accuracy, but also increasing cost, we have the simple identity | ||
| #' \deqn{E[P_L] = E[P_0] + \sum_{l=1}^L E[P_l-P_{l-1}],} | ||
| #' and therefore we can use the following unbiased estimator for \eqn{E[P_L]}, | ||
| #' \deqn{N_0^{-1} \sum_{n=1}^{N_0} P_0^{(0,n)} + \sum_{l=1}^L \{ N_l^{-1} \sum_{n=1}^{N_l} (P_l^{(l,n)} - P_{l-1}^{(l,n)}) \}} | ||
| #' with the inclusion of the level \eqn{l} in the superscript \eqn{(l,n)} indicating that the samples used at each level of correction are independent. | ||
| #' \deqn{N_0^{-1} \sum_{n=1}^{N_0} P_0^{(0,n)} + \sum_{l=1}^L \left\{ N_l^{-1} \sum_{n=1}^{N_l} \left(P_l^{(l,n)} - P_{l-1}^{(l,n)}\right) \right\}} | ||
| #' where \eqn{N_l} samples are produced at level \eqn{l}. | ||
| #' The inclusion of the level \eqn{l} in the superscript \eqn{(l,n)} indicates that the samples used at each level of correction are independent. | ||
| #' | ||
| #' Set \eqn{C_0}, and \eqn{V_0} to be the cost and variance of one sample of \eqn{P_0}, and \eqn{C_l, V_l} to be the cost and variance of one sample of \eqn{P_l - P_{l-1}}, then the overall cost and variance of the multilevel estimator is \eqn{\sum_{l=0}^L N_l C_l} and \eqn{\sum_{l=0}^L N_l^{-1} V_l}, respectively. | ||
| #' | ||
|
|
@@ -26,11 +27,11 @@ | |
| #' @author Tigran Nagapetyan <[email protected]> | ||
| #' | ||
| #' @references | ||
| #' M.B. Giles. Multilevel Monte Carlo path simulation. \emph{Operations Research}, 56(3):607-617, 2008. | ||
| #' Giles, M.B. (2008) 'Multilevel Monte Carlo Path Simulation', \emph{Operations Research}, 56(3), pp. 607–617. Available at: \url{https://doi.org/10.1287/opre.1070.0496}. | ||
| #' | ||
| #' M.B. Giles. Multilevel Monte Carlo methods. \emph{Acta Numerica}, 24:259-328, 2015. | ||
| #' Giles, M.B. (2015) 'Multilevel Monte Carlo methods', \emph{Acta Numerica}, 24, pp. 259–328. Available at: \url{https://doi.org/10.1017/S096249291500001X}. | ||
| #' | ||
| #' S. Heinrich. Monte Carlo complexity of global solution of integral equations. \emph{Journal of Complexity}, 14(2):151-175, 1998. | ||
| #' Heinrich, S. (1998) 'Monte Carlo Complexity of Global Solution of Integral Equations', \emph{Journal of Complexity}, 14(2), pp. 151–175. Available at: \url{https://doi.org/10.1006/jcom.1998.0471}. | ||
| #' | ||
| #' @param Lmin | ||
| #' the minimum level of refinement. Must be \eqn{\ge 2}. | ||
|
|
@@ -49,11 +50,12 @@ | |
| #' | ||
| #' The user supplied function should return a named list containing one element named \code{sums} and second named \code{cost}, where: | ||
| #' \describe{ | ||
| #' \item{\code{sums}}{is a vector of length two \eqn{(\sum Y_i, \sum Y_i^2)} where \eqn{Y_i} are iid simulations with expectation \eqn{E[P_0]} when \eqn{l=0} and expectation \eqn{E[P_l-P_{l-1}]} when \eqn{l>0}.} | ||
| #' \item{\code{cost}}{is a scalar with the cost of the number of paths simulated.} | ||
| #' \item{\code{sums}}{is a vector of length two \eqn{\left(\sum Y_i, \sum Y_i^2\right)} where \eqn{Y_i} are iid simulations with expectation \eqn{E[P_0]} when \eqn{l=0} and expectation \eqn{E[P_l-P_{l-1}]} when \eqn{l>0}.} | ||
| #' \item{\code{cost}}{is a scalar with the total cost of the paths simulated. | ||
| #' For example, in the financial options samplers included in this package, this is calculated as \eqn{NM^l}, where \eqn{N} is the number of paths requested in the call to the user function \code{mlmc_l}, \eqn{M} is the refinement cost factor (\eqn{M=2} for \code{\link[=mcqmc06_l]{mcqmc06_l()}} and \eqn{M=4} for \code{\link[=opre_l]{opre_l()}}), and \eqn{l} is the level being sampled.} | ||
| #' } | ||
| #' | ||
| #' See the function (and source code of) \code{\link[=opre_l]{opre_l()}} in this package for an example of a user supplied level sampler. | ||
| #' See the function (and source code of) \code{\link[=opre_l]{opre_l()}} and \code{\link[=mcqmc06_l]{mcqmc06_l()}} in this package for an example of user supplied level samplers. | ||
| #' @param alpha | ||
| #' the weak error, \eqn{O(2^{-\alpha l})}. | ||
| #' Must be \eqn{> 0} if specified. | ||
|
|
@@ -74,13 +76,13 @@ | |
| #' @return A list containing: \describe{ | ||
| #' \item{\code{P}}{The MLMC estimate;} | ||
| #' \item{\code{Nl}}{A vector of the number of samples performed on each level;} | ||
| #' \item{\code{Cl}}{Cost of samples at each level.} | ||
| #' \item{\code{Cl}}{Per sample cost at each level.} | ||
| #' } | ||
| #' | ||
| #' @examples | ||
| #' mlmc(2, 6, 1000, 0.01, opre_l, gamma=1, option=1) | ||
| #' mlmc(2, 6, 1000, 0.01, opre_l, option = 1) | ||
| #' | ||
| #' mlmc(2, 10, 1000, 0.01, mcqmc06_l, gamma=1, option=1) | ||
| #' mlmc(2, 10, 1000, 0.01, mcqmc06_l, option = 1) | ||
| #' | ||
| #' @importFrom parallel mcmapply | ||
| #' @importFrom stats lm | ||
|
|
@@ -93,15 +95,18 @@ mlmc <- function(Lmin, Lmax, N0, eps, mlmc_l, alpha = NA, beta = NA, gamma = NA, | |
| if(Lmax<Lmin) { | ||
| stop("must have Lmax >= Lmin.") | ||
| } | ||
| if(N0<=0 || eps<=0 || gamma <= 0){ | ||
| stop("N0, eps and gamma must all be greater than zero.") | ||
| if(N0<=0 || eps<=0){ | ||
| stop("N0 and eps must be greater than zero.") | ||
| } | ||
| if(!is.na(alpha) && alpha<=0) { | ||
| stop("if specified, alpha must be greater than zero. Set alpha to NA to automatically estimate.") | ||
| } | ||
| if(!is.na(beta) && beta<=0) { | ||
| stop("if specified, beta must be greater than zero. Set beta to NA to automatically estimate.") | ||
| } | ||
| if(!is.na(gamma) && gamma<=0) { | ||
| stop("if specified, gamma must be greater than zero. Set gamma to NA to automatically estimate.") | ||
| } | ||
|
|
||
| # initialise the MLMC run | ||
| est.alpha <- is.na(alpha) | ||
|
|
||
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Oops, something went wrong.