Implements Antitetich Var support

MontecarloPricing/UEOptPriceMC.m

@@ -1,11 +1,11 @@
 function [optionPrice, priceCI] = UEOptPriceMC(spotPrice, strike, rate, TTM, ...
-    putFlag, priceModel, modelParams, nSims)
-% UEOptPriceMC Calculates the price of a European option using Monte Carlo 
-% simulation.
+    putFlag, priceModel, modelParams, nSims, flagAV)
+% UEOptPriceMC Calculates the price of a European option using Monte Carlo simulation.
 %
-% This function simulates the final asset price using either the Merton or 
-% Kou jump-diffusion model and calculates the European option price and a 
-% confidence% interval for the price estimate.
+% This function simulates the final asset price using either the Merton or Kou
+% jump-diffusion model and calculates the European option price and a confidence
+% interval for the price estimate. It optionally supports the use of antithetic
+% variates to reduce variance.
 %
 % INPUTS:
 %   spotPrice   - Current price of the underlying asset.
@@ -16,22 +16,42 @@
 %   priceModel  - String specifying the asset price model ('Merton' or 'Kou').
 %   modelParams - Struct containing parameters for the chosen price model.
 %   nSims       - Number of Monte Carlo simulations to run.
+%   flagAV      - Boolean flag indicating whether to compute antithetic variates (true/false).
 %
 % OUTPUTS:
 %   optionPrice - Estimated price of the European option.
 %   priceCI     - 95% confidence interval for the option price estimate.
 %
-% The function first adds the 'PriceProcesses' directory to the MATLAB path
-% to access necessary pricing models. It then adjusts the putFlag for use 
-% in payout calculation, simulates the asset prices at maturity according 
-% to the specified model, and finally calculates the option price and 
-% confidence interval based on the discounted payoffs.
+% DESCRIPTION:
+%   The function first adds the 'PriceProcesses' directory to the MATLAB path
+%   to access necessary pricing models. It then adjusts the putFlag for use 
+%   in payout calculation, simulates the asset prices at maturity according 
+%   to the specified model, and finally calculates the option price and 
+%   confidence interval based on the discounted payoffs. The optional antithetic
+%   variates feature is used to reduce variance in Monte Carlo simulations.
+%
+% EXAMPLES:
+%   % Example with Merton model
+%   modelParams = struct('sigmaD', 0.2, 'muJ', -0.1, 'sigmaJ', 0.1, 'lambda', 0.2);
+%   [optionPrice, priceCI] = UEOptPriceMC(100, 105, 0.05, 1, false, 'Merton', modelParams, 10000, true);
+%   % This example estimates the price of a call option using 10,000 simulations
+%   % with antithetic variates enabled.
+%
+%   % Example with Kou model
+%   modelParams = struct('sigmaD', 0.2, 'lambda', 0.1, 'lambdaP', 10, 'lambdaN', 20, 'p', 0.6);
+%   [optionPrice, priceCI] = UEOptPriceMC(100, 105, 0.05, 1, true, 'Kou', modelParams, 10000, false);
+%   % This example estimates the price of a put option using 10,000 simulations
+%   % without antithetic variates.
 
 % Add the directory containing the price models to the MATLAB path
 addpath(genpath(fullfile('..', 'PriceProcesses')));
 
+if nargin < 9
+    flagAV = false;
+end
+
 % Convert the putFlag into a numerical value that will be used in the payoff calculation
-if (putFlag == true)
+if putFlag
     putFlag = -1;  % For a put option, payoff involves max(K - S, 0)
 else
     putFlag = 1;   % For a call option, payoff involves max(S - K, 0)
@@ -41,19 +61,29 @@
 switch priceModel
     case 'Merton'
         % Simulate asset prices at maturity using the Merton jump-diffusion model
-        X_T = MertonProcess(modelParams.sigmaD, modelParams.lambda, ...
-            modelParams.muJ, modelParams.sigmaJ, TTM, 1, nSims);
+        [X_T, X_T_AV] = MertonProcess(modelParams.sigmaD, modelParams.muJ, ...
+            modelParams.sigmaJ, modelParams.lambda, TTM, 1, nSims, flagAV);
     case 'Kou'
         % Simulate asset prices at maturity using the Kou jump-diffusion model
-        X_T = KouProcess(modelParams.sigmaD, modelParams.lambda, ...
+        [X_T, X_T_AV] = KouProcess(modelParams.sigmaD, modelParams.lambda, ...
                             modelParams.lambdaP, modelParams.lambdaN, ...
-                            modelParams.p, TTM, 1, nSims);
+                            modelParams.p, TTM, 1, nSims, flagAV);
 end
 
-% Extract the last simulation result as the final asset prices
-X_T = X_T(:,end);
+% Extract the simulated asset prices at maturity
+X_T = X_T(:, end);
+discountedPayoff = exp(-rate*TTM) * max(0, putFlag * (spotPrice * exp(X_T) - strike));
+
+if flagAV
+    % Compute the antithetic payoffs
+    X_T_AV = X_T_AV(:, end);
+    discountedPayoff_AV = exp(-rate*TTM) * max(0, putFlag * (spotPrice * exp(X_T_AV) - strike));
+
+    % Calculate the option price and confidence interval using both original and antithetic payoffs
+    [optionPrice, ~, priceCI] = normfit((discountedPayoff + discountedPayoff_AV) / 2);
+else
+    % Calculate the option price and confidence interval using only the original payoffs
+    [optionPrice, ~, priceCI] = normfit(discountedPayoff);
+end
 
-% Calculate the European option price and confidence interval using
-% the normal approximation of the mean of discounted payoffs
-[optionPrice, ~, priceCI] = normfit( exp(-rate*TTM) * max(0, putFlag*(spotPrice*exp(X_T) - strike)) );
 end

PriceProcesses/KouProcess.m

@@ -1,6 +1,6 @@
-function [X_t, t_i] = KouProcess(sigmaD, lambda, lambdaP, lambdaN, p, T, nTimeSteps, nProcesses)
+function [X_t, X_t_AV, t_i] = KouProcess(sigmaD, lambda, lambdaP, lambdaN, p, T, nTimeSteps, nProcesses, flagAV)
 % KOUPROCESS Simulates paths of a Kou jump-diffusion process.
-% This function generates a matrix of simulated values from a Kou jump-diffusion
+% This function generates matrices of simulated values from a Kou jump-diffusion
 % process, which is useful for modeling financial data with both continuous
 % diffusive behaviors and discontinuous jump behaviors.
 %
@@ -13,10 +13,13 @@
 %   T            - Total time span of the simulation.
 %   nTimeSteps   - Number of time steps in the simulation.
 %   nProcesses   - Number of independent processes to simulate.
+%   flagAV       - Boolean flag indicating whether to compute antithetic variates (true/false).
 %
 % OUTPUTS:
 %   X_t          - A matrix of size (nProcesses x (nTimeSteps+1)) where each row represents
-%                  a simulated path of the process.
+%                  a simulated path of the Kou model.
+%   X_t_AV       - A matrix of the same size as X_t, containing the antithetic paths.
+%                  If flagAV is false, this will be returned as an empty matrix.
 %   t_i          - A vector of length (nTimeSteps+1) representing the discrete time points
 %                  at which the process is evaluated.
 %
@@ -25,15 +28,31 @@
 %   and a jump component characterized by a double exponential distribution. Positive jumps
 %   follow one exponential distribution (with rate lambdaP), while negative jumps follow another
 %   exponential distribution (with rate lambdaN), making this model asymmetric and suitable for
-%   modeling assets with skewed jump behavior.
+%   modeling assets with skewed jump behavior. The optional antithetic variates feature is used
+%   to reduce variance in Monte Carlo simulations.
+%
+% EXAMPLES:
+%   [X_t, X_t_AV, t_i] = KouProcess(0.2, 0.1, 10, 20, 0.6, 1, 100, 1000, true);
+%   % This example generates 1000 simulated paths of a Kou process over 1 year with
+%   % 100 time steps, a diffusion volatility of 0.2, a jump intensity of 0.1,
+%   % positive and negative jump rates of 10 and 20 respectively, and returns
+%   % antithetic paths.
 
+if nargin < 9
+    flagAV = false;
+end
 
 % Calculate the time increment.
 dt = T / nTimeSteps;
 t_i = 0:dt:T;
 
-% Initialize the matrix to store the process values. First column is zero (initial value).
+% Initialize the matrices to store the process values. First column is zero (initial value).
 X_t = zeros(nProcesses, nTimeSteps+1);
+if flagAV
+    X_t_AV = zeros(nProcesses, nTimeSteps+1);
+else
+    X_t_AV = [];
+end
 
 % Define the characteristic exponent Psi of the Lévy process associated with the Kou model.
 Psi = @(u) -0.5*(sigmaD*u)^2 + 1i*u*lambda*(p/(lambdaP-1i*u) - (1-p)/(lambdaN+1i*u));
@@ -51,27 +70,44 @@
 for j = 1:nTimeSteps
     % Update each process for diffusion and drift.
     X_t(:, j+1) = X_t(:, j) + muD*dt + sigmaD * sqrt(dt) * Z(:, j);
-
+
+    if flagAV
+        X_t_AV(:, j+1) = X_t_AV(:, j) + muD*dt - sigmaD * sqrt(dt) * Z(:, j);
+    end
+
     % Loop over each process to apply the jump component.
     for i = 1<nProcesses
         % Check if there are jumps in this time step for the current process.
-        if (Ndt(i, j) > 0)
+        if Ndt(i, j) > 0
             Y = 0;
+            Y_AV = 0; % Always generated
             % Determine the number of positive and negative jumps.
             nPositiveJumps = sum(rand(Ndt(i, j), 1) < p);
             nNegativeJumps = Ndt(i, j) - nPositiveJumps;
-            
+
             % Sum the jump sizes using exponentially distributed magnitudes.
-            if (nPositiveJumps > 0)
-                Y = sum(icdf('Exponential', rand(nPositiveJumps, 1), 1/lambdaP));
+            if nPositiveJumps > 0
+                u = rand(nPositiveJumps, 1);
+                Y = sum(icdf('Exponential', u, 1/lambdaP));
+                if flagAV
+                    Y_AV = sum(icdf('Exponential', 1-u, 1/lambdaP));
+                end
             end
-            if (nNegativeJumps > 0)
-                Y = Y - sum(icdf('Exponential', rand(nNegativeJumps, 1), 1/lambdaN));
+            if nNegativeJumps > 0
+                u = rand(nNegativeJumps, 1);
+                Y = Y - sum(icdf('Exponential', u, 1/lambdaN));
+                if flagAV
+                    Y_AV = Y_AV - sum(icdf('Exponential', 1-u, 1/lambdaN));
+                end
             end
-            
+
             % Update the process value by the total jump magnitude.
             X_t(i, j+1) = X_t(i, j+1) + Y;
+            if flagAV
+                X_t_AV(i, j+1) = X_t_AV(i, j+1) + Y_AV;
+            end
         end
     end
 end
+
 end

PriceProcesses/MertonProcess.m

@@ -1,8 +1,8 @@
-function [X_t, t_i] = MertonProcess(sigmaD, muJ, sigmaJ, lambda, T, nTimeSteps, nProcesses)
+function [X_t, X_t_AV, t_i] = MertonProcess(sigmaD, muJ, sigmaJ, lambda, T, nTimeSteps, nProcesses, flagAV)
 % MERTONPROCESS Simulates paths of a Merton jump-diffusion process.
-% This function generates a matrix of simulated values from a Merton jump-diffusion
+% This function generates matrices of simulated values from a Merton jump-diffusion
 % process, which incorporates both Gaussian diffusion elements and a jump process
-% characterized by normally distributed jumps.
+% characterized by normally distributed jumps. It optionally supports the use of antithetic variates.
 %
 % INPUTS:
 %   sigmaD       - Volatility of the diffusion part of the stock price.
@@ -12,18 +12,26 @@
 %   T            - Total time span of the simulation.
 %   nTimeSteps   - Number of time steps in the simulation.
 %   nProcesses   - Number of independent process paths to simulate.
+%   flagAV       - Boolean flag indicating whether to compute antithetic variates (true/false).
 %
 % OUTPUTS:
 %   X_t          - A matrix of size (nProcesses x (nTimeSteps+1)) where each row represents
 %                  a simulated path of the Merton model.
+%   X_t_AV       - A matrix of the same size as X_t, containing the antithetic paths.
+%                  If flagAV is false, this will be returned as an empty matrix.
 %   t_i          - A vector of length (nTimeSteps+1) representing the discrete time points
 %                  at which the process is evaluated.
 %
 % DESCRIPTION:
 %   The Merton jump-diffusion model includes both a continuous Gaussian diffusion part
 %   and a jump component characterized by a Poisson process with normally distributed jump sizes.
 %   This allows modeling assets with occasional large jumps in price in addition to the standard
-%   continuous diffusive behavior.
+%   continuous diffusive behavior. The optional antithetic variates feature is used to reduce variance
+%   in Monte Carlo simulations.
+
+if nargin < 8
+    flagAV = false;
+end
 
 % Calculate the time increment.
 dt = T / nTimeSteps;
@@ -49,13 +57,21 @@
 
 % Initialize a matrix to accumulate jump sizes.
 Y = zeros(nProcesses, nTimeSteps);
+if flagAV
+    Y_AV = zeros(nProcesses, nTimeSteps);
+end
 
 % Loop over each process and time step to accumulate jumps.
 for i = 1:nProcesses
     for j = 1:nTimeSteps
         if Ndt(i, j) > 0
             % Sum up the jump sizes, where jump magnitudes are normally distributed.
             Y(i, j) = sum(muJ + sigmaJ * randn(1, Ndt(i, j)));
+
+            if flagAV
+                % Use antithetic variates for the jump sizes.
+                Y_AV(i, j) = sum(muJ - sigmaJ * randn(1, Ndt(i, j)));
+            end
         end
     end
 end
@@ -64,9 +80,22 @@
 X_t = muW * t_i .* ones(nProcesses, nTimeSteps) ...
         + B_t ...
         + cumsum(Y, 2);
-        
+
 % Prepend zeros to represent the starting value of the process.
 X_t = cat(2, zeros(nProcesses, 1), X_t);
 t_i = cat(2, 0, t_i);
 
+if flagAV
+    % Compute the antithetic paths with the opposite sign of the diffusion component.
+    X_t_AV = muW * t_i .* ones(nProcesses, nTimeSteps) ...
+        - B_t ...
+        + cumsum(Y_AV, 2);
+
+    % Prepend zeros to represent the starting value of the antithetic process.
+    X_t_AV = cat(2, zeros(nProcesses, 1), X_t_AV);
+else
+    % Return an empty matrix if antithetic variates are not requested.
+    X_t_AV = [];
+end
+
 end