clear all;

global tau;

%----------------------------- parameters
tau = 1; %accrual time
alpha=5; %start peg of swap
beta=9; %end peg of swap
dt=1; %time-step of monte carlo simulation
N=10; %scenarios
L=10000;
%----------------------------------------------

function ret = volatility(Tj,T_o)
	a = -0.05;
	b = 0.5;
	c = 1.5;
	d = 0.15;
	kj = 1;
	ret=kj * ((a + b * (Tj - T_o)) * exp(-c * (Tj - T_o)) + d);
endfunction

function ret = corr(Tj,Tk)
	beta = 0.1;
	ret=exp(-beta * abs(Tj - Tk));
endfunction

%based on eq. 6.32 
function ret = GetSwapRate(F,alpha,beta)
	global tau;
	tmp_sum=0;
	SR=1;
	tmp=1;
	for j=alpha+1:beta,
		tmp=tmp*(1/(1+tau*F(j-1)));
	end
	SR=1-tmp; %numerator 
	for i=alpha+1:beta,
		tmp=1;
		for j=alpha+1:i,
			tmp=tmp*(1/(1+tau*F(j-1)));	  
		end
		tmp_sum=tmp_sum + (tau*tmp);
	end
	SR=SR/tmp_sum;
	ret=SR;
endfunction

%load -force -ascii p_matrix.txt 
%discount curve as seen today
P=[1.000000
0.966736
0.934579
0.903492
0.873439
0.844385
0.816298
0.789145
0.762895
0.737519
0.712986
0.689270
0.666342
0.644177
0.622750
0.602035
0.582009
0.562649
0.543934
0.525841
0.508349
0.491440
0.475093
0.459290
0.444012
0.429243
0.414964
0.401161
0.387817
0.374917
0.362446
0.350390
0.338735
0.327467
0.316574
0.306044
0.295864
0.286022
0.276508
0.267311
0.258419
0.249823]; 

F=zeros(size(P,1)-1,1);
T=0:tau:20.5;
for i=1:size(F),
	F(i)=(P(i)/P(i+1)-1)/tau;
end

%F=F*0+0.05;

SR=GetSwapRate(F,alpha,beta);
SR

nF=size(F); %no. of froward rates
rootdT=dt^0.5;
%nFswap=beta-alpha;
nFswap=beta-1;
%form a correlation matrix rho which represents 
%the correlation between forward rates
rho=zeros(nFswap,nFswap);
tmp_row=0;
tmp_col=0;
%for i=alpha:beta-1,
for i=1:beta-1,
	tmp_row=tmp_row+1;
	tmp_col=0;
	%for j=alpha:beta-1,
	for j=1:beta-1,
		tmp_col=tmp_col+1;
		rho(tmp_row,tmp_col)=corr(T(i),T(j));
	end
end
chol_rho=chol(rho);
chol_rho
mu=zeros(nFswap,1);
K=SR; %strike of swaption is ATM swap rate

payoff_sum=0;
payoff_sum2=0;
flag_antithetic=-1;
for scenario=1:N, 
	logF=log(F); 
	cases=T(alpha)/dt; %only the forward rates within swap need to be simulated
	if flag_antithetic == -1, %use antithetic MC
		tmprndvec=randn(cases+10,size(rho,1));
		%tmprndvec(:,2:size(rho,1))=repmat(tmprndvec(:,1),1,size(rho,1)-1);
		dZ=tmprndvec*chol_rho;
		%dZ=tmprndvec;
		dZ=dZ*rootdT;
	else
		dZ=-dZ;
	endif
	flag_antithetic=flag_antithetic*-1;
	tmpcount=1;
	t=0; %current time
	while t<=T(alpha), %simulate the path for each forward rate, till beginning of swap
		%t=t+dt;
		t;
		tmp_col=0;
		startk=floor(t/tau)+1;
		%for k=alpha:beta-1,
		for k=startk:beta-1,
			drift_sum=0;
			%for j=alpha:k,
			for j=startk:k+1,
				%tmp=(corr(T(k),T(j))*tau*volatility(T(j),t)*exp(logF(j)) );
				%tmp=tmp/(1+tau*exp(logF(j)));
				
				tmp=(corr(T(k),T(j))*tau*volatility(T(j),t)*F(j) );
				%tmp=(tau*0.2*F(j) );
				tmp=tmp/(1+tau*F(j));
				
				drift_sum = drift_sum+tmp;
			end
			tmp_col=tmp_col+1;
			%logF(k)=logF(k)+ volatility(T(k),t)*drift_sum*dt - (volatility(T(k),t)^2*dt)/2 + volatility(T(k),t)*dZ(tmpcount,tmp_col);
			%logF(k+1)=logF(k+1)+ 0.2*drift_sum*dt - (0.2^2*dt)/2 + 0.2*dZ(tmpcount,tmp_col);
			sigma_Tk_t=volatility(T(k),t);
			logF(k+1)=logF(k+1)+ sigma_Tk_t*drift_sum*dt - (sigma_Tk_t*sigma_Tk_t*dt)/2 + sigma_Tk_t*dZ(tmpcount,tmp_col);
		end
		tmpcount=tmpcount+1;
		t=t+dt;
	endwhile
	Ffinal=exp(logF); %after simulating path of each forward rate, final F vector in alpha forward measure
	SR=GetSwapRate(Ffinal,alpha,beta);
	tmp_sum=0;
	for i=alpha+1:beta,
		tmp=1;
		for j=alpha+1:i,
			tmp=tmp*(1/(1+tau*Ffinal(j-1)));
		end
		tmp_sum=tmp_sum+tau*tmp;
	end
	df=1;
	for i=1:alpha-1,
		df=df*(1/(1+tau*Ffinal(j)));
	end
	%payoff=P(alpha)*max(SR-K,0)*tmp_sum;
	payoff=df*max(SR-K,0)*tmp_sum;
	payoff_sum=payoff_sum+payoff;
	
	%payoff2=L*df*tau*max(Ffinal(alpha)-0.05,0);
	%payoff_sum2=payoff_sum2+payoff2;
end %for each .. scenario

payoff=payoff_sum/N;
%payoff2=payoff_sum2/N
swaption_price=payoff*100;
tmpstr=sprintf('swap start=%f   \nswap end=%f      \nswaption price=%f',T(alpha),T(beta),swaption_price);
disp(tmpstr);