% uppg5_13.m % % An example of variance reduction % % Goal : compute % % pi_t = exp(-r*(T-t))*E[max(Z(T)-K,0)] % % where % % Z(t) = (1/d)*sum(S_i(t)) % % is the average of d stock prices which in the risk neutral % formulation follows the SDE:s % % dS_i(t) = r*S_i(t) dt + sigma_i*S_i(t) dW_i % % The volatilities are % % sigma_i = 0.2*(2+sin(i)), i = 1,...,d % % and the Weiner processes are correlated as % % E[dW_i(t)*dW_j(t)] = exp(-2|i-j|/d)dt % Set the state of the generator of normally distributed random % variables. This gaurantees that the sequence of pseudo random % numbers is identical in each run. randn('state',100); % ---------- Monte-Carlo parameters ---------------------------- Mv = 2.^(6:23); % Number of samples (for convergence study) Cconf = 2; % Confidence interval parameter TOL = 2e-3; % Tolerance % ---------- Problem parameters -------------------------------- disp('Problem parameters : ') d = 10 % Number of stocks r = 0.04 % Risk free interest rate T = 0.5 % Final time K = 40 % v = 0.2*(2+sin(1:d))' % Volatilities S0 = K*ones(d,1); % Initial stock prices % Correlation matrix b = exp(-2*abs((1-d):(d-1))/d)'; B = zeros(d,d); for k=1:d, B(:,k) = b((d-k+1):(2*d-k)); end; % $$$ % Try highly correlated stocks % $$$ B = 0.98*ones(d,d); B=B+0.02*spdiags(ones(d,1),0,d,d); % Display the correlation matrix disp('Correlation matrix for the noise in the stock price processes') for k=1:d, str = ''; for l=1:d, str = [str, num2str(B(k,l),'%1.2f'),' ']; end; disp(str); end; % ------ end Problem parameters -------------------------------- % Transforms d independent Wienerprocesses to d processes with % correlation matrix B, using a linear mapping. [V,L] = eig(B); trans = V*sqrt(L); % Diffusion matrix of d-dimensional Ito-SDE with independent Wiener % processes sigma = spdiags(v,0,d,d)*trans; % ---------- Control variate Z_hat setup ----------------------- Z_hat0 = (prod(S0))^(1/d); sig_hat = mean(sigma); % Calculate the drift of the 1-dim. 'approximating' SDE alpha = 0; for i=1:d, for j=1:d, for l=1:d, alpha = alpha + sigma(i,l)*sigma(j,l); end; end; end; alpha = alpha/(2*d^2); alpha = alpha + r; tmp = 0; for i=1:d, for l=1:d, tmp = tmp + sigma(i,l)^2; end; end; tmp = tmp/(2*d); alpha = alpha - tmp; % Calculate the diffusion of the 1-dim. 'approximating' SDE beta = 0; for l=1:d, beta = beta + (sum(sigma(:,l)))^2; end; beta = (1/d)*sqrt(beta); % Use Black-Scholes formula to evaluate % exp(-alpha*T)*E(max(Z_hat(T)-K,0)) g = BS_form(0,S0(1),T,K,alpha,beta) % ------ end Control variate Z_hat setup ---------------------- % ---------- Main Monte Carlo loop ---------------------------- Pi_t = []; Pi_t_anti = []; Pi_t_cv = []; Pi_t_cv_anti = []; errest1 = []; errest2 = []; errest3 = []; errest4 = []; Time1 = []; Time2 = []; Time3 = []; Time4 = []; Time_randn = []; % $$$ % Uncomment to stop at a specified relative tolerance % $$$ kk = 0 % $$$ while ((kk==0) | (errest4(end)>TOL*Pi_t_cv_anti(end))), % $$$ kk = kk+1; M = Mv(kk); for M = Mv, tic; % Generate M samples of d-dim. Wiener process at time T W_T = sqrt(T)*randn(d,M); t_randn = toc; Time_randn = [Time_randn t_randn]; tic; Z_T = zeros(1,M); Z_anti = zeros(1,M); % Use explicit solution formula for SDE to evaluate Z(T) % Z(T) is arithmetic mean of stock prices for k=1:d, tmp1 = sigma(k,:)*W_T; tmp2 = 0.5*sigma(k,:)*(sigma(k,:))'*T; Z_T = Z_T + S0(k)*exp(r*T + tmp1 - tmp2); %antithetic variables Z_anti = Z_anti + S0(k)*exp(r*T - tmp1 - tmp2); end; Z_T = Z_T/d; Z_anti = Z_anti/d; gval = exp(-r*T)*max(Z_T-K,0); gval_anti = exp(-r*T)*max(Z_anti-K,0); t1 = toc; % Method 1 : Monte Carlo without variance reduction tic; errest1 = [errest1 Cconf*std(gval)/sqrt(M)]; %check = errest1(end) - Cconf*sqrt(sum((gval-mean(gval)).^2)/(M-1))/sqrt(M) pi_t = mean(gval); % Sample average t2 = toc; % Method 2 : Monte Carlo using antithetic variables tic; gval2 = mean([gval;gval_anti]); errest2 = [errest2 Cconf*std(gval2)/sqrt(M)]; pi_t_anti = mean(gval2); t3 = toc; % ------ Control variate computations -------------------- tic; tmp1 = sig_hat*W_T; tmp2 = 0.5*sig_hat*sig_hat'*T; Z_hat = Z_hat0*exp(alpha*T + tmp1 - tmp2); Z_hat_anti = Z_hat0*exp(alpha*T - tmp1 - tmp2); gval3 = gval - exp(-alpha*T)*max(Z_hat-K,0); gval3_anti = gval_anti - exp(-alpha*T)*max(Z_hat_anti-K,0); t4 = toc; % Method 3 : Monte Carlo using control variate tic; errest3 = [errest3 Cconf*std(gval3)/sqrt(M)]; pi_t_cv = g + mean(gval3); t5 = toc; % Method 4 : Monte Carlo using control variate and antithetic % variables tic; gval4 = mean([gval3; gval3_anti]); errest4 = [errest4 Cconf*std(gval4)/sqrt(M)]; pi_t_cv_anti = g + mean(gval4); t6 = toc; % Output and storage disp(['pi_t = ',num2str(pi_t),... ' , pi_t_anti = ',num2str(pi_t_anti),... ' , pi_t_cv = ',num2str(pi_t_cv),... ' , pi_t_cv_anti = ',num2str(pi_t_cv_anti)]); Pi_t = [Pi_t pi_t]; Time1 = [Time1 (t_randn+0.5*t1+t2)]; Pi_t_anti = [Pi_t_anti pi_t_anti]; Time2 = [Time2 (t_randn+t1+t3)]; Pi_t_cv = [Pi_t_cv pi_t_cv]; Time3 = [Time3 (t_randn+0.5*t1+t2+0.5*t4+t5)]; Pi_t_cv_anti = [Pi_t_cv_anti pi_t_cv_anti]; Time4 = [Time4 (t_randn+t1+t3+t4+t6)]; end; % $$$ % Uncomment to stop at a specified relative tolerance % $$$ Mv = Mv(1:kk); % ---------- Plots below ----------------------------------- figure, p=plot(log2(Mv),Pi_t,'k-o',log2(Mv),Pi_t_anti,'r-x',... log2(Mv),Pi_t_cv,'b-d',log2(Mv),Pi_t_cv_anti,'m-*'); legend(p,'standard','antithetic','control variates',... 'control variates and antithetic'); xlabel('log_2( Nr. of Samples )','fontsize',14); ylabel('\pi_t','fontsize',14,'rotation',0); figure, ref_sol = Pi_t_cv_anti(end); p=plot(log2(Mv),log2(abs(Pi_t-ref_sol)),'k-o',... log2(Mv),log2(abs(Pi_t_anti-ref_sol)),'r-x',... log2(Mv),log2(abs(Pi_t_cv-ref_sol)),'b-d',... log2(Mv(1:end-2)),log2(abs(Pi_t_cv_anti(1:end-2)-ref_sol)),'m-*'); legend(p,'standard','antithetic','control variates',... 'control variates and antithetic'); xlabel('log_2( Nr. of Samples )','fontsize',14); ylabel('log_2(|\pi_t-ref.val|)','fontsize',14,'rotation',0); axis image; grid on; figure, p=plot(log2(Mv),log2(errest1),'k-o',... log2(Mv),log2(errest2),'r-x',... log2(Mv),log2(errest3),'b-d',... log2(Mv),log2(errest4),'m-*'); legend(p,'standard','antithetic','control variates',... 'control variates and antithetic'); xlabel('log_2( Nr. of Samples )','fontsize',14); ylabel('log_2(|errest|)','fontsize',14,'rotation',0); axis image; grid on; Time1 = cumsum(Time1); Time2 = cumsum(Time2); Time3 = cumsum(Time3); Time4 = cumsum(Time4); T1 = Time1(min(find(errest1