%% EXERCISE 1 - STYLISED FACTS

% Pavol Jurca, 2013

% The aim of this exercise is to show that standard assumptions 
% of normality, stationarity and simple correlation structure 
% are not fulfilled in real world

%% Data preparation

clearvars
clc

% Load historical stock prices for S&P 500 and FTSE 100
% Time series from 1 January 1985 to 31 August 2013

prices=xlsread('stocks.xls');

% calculate relative (continuously compounded) returns
returns=diff(log(prices));
returns1=returns(:,1);
returns2=returns(:,2);

figure(1)
plot(returns1)

%% Non-normality

% Plot histogram with 100 bins together with normal fit
histfit(returns1,100)

% Jarque-Bera test of normality 
% (returns 1 if the null hypothesis of normality is rejected)
[h,p]=jbtest(returns1)

% Quantile-quantile plot of the sample quantiles returns1 
% versus theoretical quantiles from a normal distribution
figure(2)
subplot(1,2,1)
qqplot(returns1)

% Open the distribution fitting application
%dfittool(returns1)

% QQ plot versus t distribution with d.o.f.=2.75
subplot(1,2,2)
qqplot(returns1,trnd(2.88,length(returns1),1))


%% Fat tails

figure(3)
[F,yi] = ecdf(returns1);   % empirical CDF
stairs(yi,F,'r'); hold on; 
axis([-0.11 -0.015 0 0.1]);

plot(yi,normcdf(yi,mean(returns1),std(returns1)),'-'); % fitted normal
hold off;

legend('Fitted Normal CDF','Empirical CDF','location','northwest');

%% Volatility clustering

% Autocorrelation function of ordinary as well as squared residuals
figure(4)
subplot(1,2,1)
autocorr(returns1)
subplot(1,2,2)
autocorr(returns1.^2)

% Standard deviation in one-year periods
years=1985:1:2012;
returns1_by_years=reshape(returns1(1:7000,:),250,28);
returns2_by_years=reshape(returns2(1:7000,:),250,28);

figure(5)
bar(years,std(returns1_by_years))


%% Nonlinear dependence

% Multivariate distribution

% Comparison to multivariate normal and multivariate Student t
sim_normal=mvnrnd(mean(returns),cov(returns),length(returns1));
sim_t=mvtrnd(corr(returns),2.88,length(returns1));
sim_t1=sim_t(:,1).*(std(returns1)*sqrt(0.88/2.88));
sim_t2=sim_t(:,2).*(std(returns2)*sqrt(0.88/2.88));

figure(6)

subplot(1,3,1)
scatter(returns1, returns2);
axis([-0.2 0.2 -0.2 0.2]);

subplot(1,3,2)
scatter(sim_normal(:,1),sim_normal(:,2));
axis([-0.2 0.2 -0.2 0.2]);

subplot(1,3,3)
scatter(sim_t1,sim_t2);
axis([-0.2 0.2 -0.2 0.2]);

% Time varying correlations

for i=1:28
    correlation_by_year(i)=corr(returns1_by_years(:,i),returns2_by_years(:,i));
end    

figure(7)
bar(years,correlation_by_year)

% Exceedance correlations

% Calculation of exceedance correlations for original data
p=0.03:0.01:0.97;
exceedance_correlation=ones(size(p));
quantile1=quantile(returns1,p);
quantile2=quantile(returns2,p);

for i=1:length(p)
   if p(i)<=0.5
        indexes=and(returns1<=quantile1(i),returns2<=quantile2(i));
    else     
        indexes=and(returns1>quantile1(i),returns2>quantile2(i));
    end    
    returns1_censored=returns1(indexes==1);
    returns2_censored=returns2(indexes==1);
    exceedance_correlation(i)=corr(returns1_censored,returns2_censored);
end

figure(8)
plot(p,exceedance_correlation,'-r')
hold on

% Calculation of exceedance correlations for multivariate normal data
simulated_normal=mvnrnd(mean(returns),cov(returns),200000);
returns1_sim=simulated_normal(:,1);
returns2_sim=simulated_normal(:,2);

quantile1=quantile(returns1_sim,p);
quantile2=quantile(returns2_sim,p);

for i=1:length(p)
    if p(i)<=0.5
        indexes=and(returns1_sim<=quantile1(i),returns2_sim<=quantile2(i));
    else     
        indexes=and(returns1_sim>quantile1(i),returns2_sim>quantile2(i));
    end    
    returns1_censored=returns1_sim(indexes==1);
    returns2_censored=returns2_sim(indexes==1);
    exceedance_correlation_sim(i)=corr(returns1_censored,returns2_censored);
end

plot(p,exceedance_correlation_sim,'-.g')
hold on

% Calculation of exceedance correlations for multivariate Student t data
simulated_t=mvtrnd(cov(returns),3,1000000);
returns1_t=simulated_t(:,1);
returns2_t=simulated_t(:,2);

quantile1=quantile(returns1_t,p);
quantile2=quantile(returns2_t,p);

for i=1:length(p)
    if p(i)<=0.5
        indexes=and(returns1_t<=quantile1(i),returns2_t<=quantile2(i));
    else     
        indexes=and(returns1_t>quantile1(i),returns2_t>quantile2(i));
    end    
    returns1_censored=returns1_t(indexes==1);
    returns2_censored=returns2_t(indexes==1);
    exceedance_correlation_t(i)=corr(returns1_censored,returns2_censored);
end

plot(p,exceedance_correlation_t, '--b')