function [a,b]=h1(h,M) 

% option on average price of two assets obeying coupled GBM
% want sensitivity of option price with respect to one underlying 

% will find derivative using f'(x) = ( f(x+h) - f(x-h) ) / (2h) 
% using M simulations for each of f(x+h) and f(x-h) 

% all option parameters 
S10 = 15;
S20 = 30; 
r = 0.05;
sigma1 = 0.1 ;
sigma2 = 0.08; 
rho = 0.4; 
T = 1;
K = 24; 

% "f(x+h)" 
Z = randn(M,2);
V2 = rho * Z(:,1) +  sqrt(1-rho^2) * Z(:,2) ; 
S1 = (S10 + h) * exp( (r-sigma1^2/2)*T + sigma1 * sqrt(T) * Z(:,1) );
S2 = S20 * exp( (r-sigma2^2/2)*T + sigma2 * sqrt(T) * V2 ); 
A = (S1+S2)/2 ;
P1 = exp(-r*T)*(A-K).*(A>K); 

% "f(x-h)" 
Z = randn(M,2);
V2 = rho * Z(:,1) +  sqrt(1-rho^2) * Z(:,2) ; 
S1 = (S10 - h) * exp( (r-sigma1^2/2)*T + sigma1 * sqrt(T) * Z(:,1) );
S2 = S20 * exp( (r-sigma2^2/2)*T + sigma2 * sqrt(T) * V2 ); 
A = (S1+S2)/2 ;
P2 = exp(-r*T)*(A-K).*(A>K); 

delta = (P1-P2)/(2*h) ;
a =  mean(delta);
b =  std(delta)/sqrt(M);
