function VaR=MC(NbTraj,p,alpha)
NbTraj= 50000;
p= 0.5;
alpha= 0.96;
fid=fopen('SpotRates.txt','rt');
Data=fscanf(fid, '%f %f %f %f %f %f %f %f %f');
fclose(fid);
M1=Data(1:9:size(Data));
M3=Data(2:9:size(Data));
M6=Data(3:9:size(Data));
Y1=Data(4:9:size(Data));
Y2=Data(5:9:size(Data));
Y3=Data(6:9:size(Data));
Y5=Data(7:9:size(Data));
Y10=Data(8:9:size(Data));
Y30=Data(9:9:size(Data));
Data1=[M1,M3,M6,Y1,Y2,Y3,Y5,Y10,Y30];
ActualRates=Data1(size(Data1,1),:);
for i=1:9
Average(i)=mean(Data1(:,i));
Deviation(i)=std(Data1(:,i));
Data1(:,i)=(Data1(:,i)-Average(i))/Deviation(i);
end
Deviation=Deviation/sqrt(3);
[V,E]=eigs(cov(Data1),3);
PC1=V(:,1);
PC2=V(:,2);
PC3=V(:,3);
plot((1:1:9),sqrt(E(1,1))*PC1,'r',(1:1:9),sqrt(E(2,2))*PC2,'b',...
(1:1:9),sqrt(E(3,3))*PC3,'m');
Weight1=p;
Weight2=1-Weight1;
CouponRate1=0.06;
FaceValue1=1000;
CouponRate2=0.08;
FaceValue2=1000;
coupon1=CouponRate1*FaceValue1;
coupon2=CouponRate2*FaceValue2;
u1=zeros(NbTraj,1);
u2=zeros(NbTraj,1);
for l=1:NbTraj
u1(l)=norminv(VanDerCorput(l,3));
u2(l)=norminv(VanDerCorput(l,5));
end
Rates=zeros(9,NbTraj);
for j=1:NbTraj
for i=1:9
Rates(i,j)=ActualRates(1,i)+...
Deviation(i)*(sqrt(E(1,1))*u1(j,1)*V(i,1)*sign(V(1,1))+...
sqrt(E(2,2))*u2(j,1)*V(i,2)*sign(V(1,2)));
end
end
M6=ActualRates(1,3);
Y1=ActualRates(1,4) ;
Y2=ActualRates(1,5);
Y3=ActualRates(1,6);
Y5=ActualRates(1,7);
Y1_5=(Y1+Y2)/2;
Y2_5=(Y2+Y3)/2;
Y3_5=(3*Y3+Y5)/4;
Y4=(Y3+Y5)/2;
Y4_5=(3*Y5+Y3)/4;
price1=(coupon1/2)/(1+(M6/100))^(1/2)+(coupon1/2)/(1+(Y1/100))^1+...
(coupon1/2)/(1+(Y1_5/100))^(3/2)+(coupon1/2)/(1+(Y2/100))^2+...
(FaceValue1+(coupon1/2))/(1+(Y2_5/100))^(5/2);
price2=(coupon2/2)/(1+(M6/100))^(1/2)+(coupon2/2)/(1+(Y1/100))^1+...
(coupon2/2)/(1+(Y1_5/100))^(3/2)+(coupon2/2)/(1+(Y2/100))^2+...
(coupon2/2)/(1+(Y2_5/100))^(5/2)+(coupon2/2)/(1+(Y3/100))^3+...
(coupon2/2)/(1+(Y3_5/100))^(7/2)+(coupon2/2)/(1+(Y4/100))^4+...
(coupon2/2)/(1+(Y4_5/100))^(9/2)+...
(FaceValue2+(coupon2/2))/(1+(Y5/100))^5;
InitialPrice=Weight1*price1+Weight2*price2;
for k=1:NbTraj
M6=Rates(3,k);
Y1=Rates(4,k);
Y2=Rates(5,k);
Y3=Rates(6,k);
Y5=Rates(7,k);
Y1_5=(Y1+Y2)/2;
Y2_5=(Y2+Y3)/2;
Y3_5=(3*Y3+Y5)/4;
Y4=(Y3+Y5)/2;
Y4_5=(3*Y5+Y3)/4;
price1n=(coupon1/2)/(1+(M6/100))^(1/2)+(coupon1/2)/(1+(Y1/100))^1+...
(coupon1/2)/(1+(Y1_5/100))^(3/2)+(coupon1/2)/(1+(Y2/100))^2+...
(FaceValue1+(coupon1/2))/(1+(Y2_5/100))^(5/2);
price2n=(coupon2/2)/(1+(M6/100))^(1/2)+(coupon2/2)/(1+(Y1/100))^1+...
(coupon2/2)/(1+(Y1_5/100))^(3/2)+(coupon2/2)/(1+(Y2/100))^2+...
(coupon2/2)/(1+(Y2_5/100))^(5/2)+(coupon2/2)/(1+(Y3/100))^3+...
(coupon2/2)/(1+(Y3_5/100))^(7/2)+(coupon2/2)/(1+(Y4/100))^4+...
(coupon2/2)/(1+(Y4_5/100))^(9/2)+...
(FaceValue2+(coupon2/2))/(1+(Y5/100))^5;
PortfolioPrice=Weight1*price1n+Weight2*price2n;
vectPrice(k,1)= PortfolioPrice;
end
vectPrice=sort(vectPrice);
VaR=(InitialPrice-vectPrice(floor(alpha/100*NbTraj)));
vectprice1(k,1)=price1n;
end
function rep=VanDerCorput(n,b)
n=100;
b=2;
bn=0;
j=0;
while n~=0
bn=bn+mod(n,b)/b^(j+1);
n=floor(n/b);
j=j+1;
end
rep=bn;
end