
%(c) Pavol Jurca, 2014

clc
Data=xlsread('fx_rates.xls'); %Data for USD	JPY	GBP	CHF	AUD	CAD, Jan 1999- Dec 2012
%plot(Data)
Returns=log(Data(2:end,:)./Data(1:(end-1),:)); %calculating relative changes

%MIXTURE OF 2 NORMAL DISTRIBUTION, 1 MAIN RISK FACTOR
obj=gmdistribution.fit(Returns(:,1),2);
Proportions = obj.PComponents;
Means=obj.mu;
Volatilities=sqrt(obj.Sigma);

if Volatilities(:,:,2)>Volatilities(:,:,1)
   hectic=2;
else 
   hectic=1;
end

%Calculating parameters for the primary risk factor
disp('Calculating parameters for the primary risk factor') 
quiet=3-hectic;
omega=Proportions(hectic)
mu2=Means(hectic)
mu1=Means(quiet)
sigma2=Volatilities(hectic)
sigma1=Volatilities(quiet)

pause




%Calculating probablities that an observation is from hectic period
h=(omega.*pdf('Normal',Returns(:,1),mu2,sigma2))./((1-omega).*pdf('Normal',Returns(:,1),mu1,sigma1)+omega.*pdf('Normal',Returns(:,1),mu2,sigma2));
scatter(Returns(:,1),h)
pause

%Calculating parameters for the secondary risk factor            
disp('Calculating parameters for peripheral risk factor') 
scenarios=[0.3 0 0 0 0 0];
for i=2:6
   N=size(Returns(:,i));
   mu_y2=sum(h.*Returns(:,i))/sum(h)
   mu_y1=sum((1-h).*Returns(:,i))/sum(1-h)
   sigma_y2=sqrt(sum(h.*(Returns(:,i)-mu_y2*ones(N)).^2)/sum(h))
   sigma_y1=sqrt(sum((1-h).*(Returns(:,i)-mu_y1*ones(N)).^2)/sum(1-h))
   rho_y2=sum(h.*(Returns(:,1)-mu2*ones(N)).*(Returns(:,i)-mu_y2*ones(N)))/(sigma2*sigma_y2*sum(h))
   rho_y1=sum((1-h).*(Returns(:,1)-mu1*ones(N)).*(Returns(:,i)-mu_y1*ones(N)))/(sigma1*sigma_y1*sum(1-h))
   scenarios(i)=rho_y2*(scenarios(1)-mu2)*sigma_y2/sigma2+mu_y2
   pause
end



pause

%MIXTURE OF 2 NORMAL DISTRIBUTION, 2 EQUIVALENT RISK FACTORS
disp('MIXTURE OF 2 NORMAL DISTRIBUTION, 2 EQUIVALENT RISK FACTORS')
obj2=gmdistribution.fit(Returns(:,1:2),2);
Cov=obj2.Sigma;
Cov1=Cov(:,:,1);
Cov2=Cov(:,:,2);
rho_1=Cov1(1,2)/sqrt(Cov1(1,1)*Cov1(2,2))
rho_2=Cov2(1,2)/sqrt(Cov2(1,1)*Cov2(2,2))

pause

%TEST OF THE SPECIFICATION OF THE MODEL FOR 1 MAIN RISK FACTOR
disp('TEST OF THE SPECIFICATION OF THE MODEL FOR 1 MAIN RISK FACTOR')
temp_obj=gmdistribution.fit(Returns(:,1),2);
LogL2=temp_obj.NlogL;
temp_obj=gmdistribution.fit(Returns(:,1),1);
LogL1=temp_obj.NlogL;
Likelihood_ratio=-2*(LogL2-LogL1)
ch2_quantil=chi2inv(0.95,3)

pause

%MIXTURE OF 3 NORMAL DISTRIBUTION, 1 RISK FACTOR
disp('MIXTURE OF 3 NORMAL DISTRIBUTION, 1 RISK FACTOR')
options = statset('MaxIter',1000);
obj3=gmdistribution.fit(Returns(:,1:1),3,'Options',options);
Sigma=sqrt(obj3.Sigma)
Proportions=obj3.PComponents
