%% Import my daily prices on 13 forex pairs (first column are dates) - period of 15 years

data = xlsread('Daily Data FX.xlsx');
dates = data(:,1);


%% Get saturdays out (Weird that is needs #1 but this seems to work)

dates = dates + datenum('31DEC1899');
saturdays = weekday(dates) == 1;
data(saturdays,:) = [];
dates(saturdays) = [];
clearvars saturdays;

%% Get prices from data vector

prices = data(:,2:end);


%% Calculate moving averages

MA10 = tsmovavg(prices,'e',10,1);
MA50 = tsmovavg(prices,'e',50,1);
MA200 = tsmovavg(prices,'e',200,1);


%% Get rid of first row

data(1,:) = [];
dates(1,:) = [];
MA10(1,:) = [];
MA50(1,:) = [];
MA200(1,:) = [];
prices(1,:) = [];

%% Get rid of last 4 rows

data(end-4:end,:) = [];
dates(end-4:end,:) = [];
MA10(end-4:end,:) = [];
MA50(end-4:end,:) = [];
MA200(end-4:end,:) = [];
prices(end-4:end,:) = [];

%% Calculate continuous returns

returns = price2ret(prices,'','cont');

%% Moments

Moments_per_pair(1,:) = mean(returns);
Moments_per_pair(2,:) = median(returns);
Moments_per_pair(3,:) = max(returns);
Moments_per_pair(4,:) = min(returns);
Moments_per_pair(5,:) = skewness(returns);
Moments_per_pair(6,:) = kurtosis(returns);

%% Logical for trending above all MA

trending_up = (prices > MA10) & (prices > MA50) & (prices > MA200);
trending_down = (prices < MA10) & (prices < MA50) & (prices < MA200);

%% Calculate returns for all observations in 1 column

returns_all_one_column = reshape(returns,[],1);

%% Calculate returns when trending conditions are met

returns(end+1,:) = 0;
returns_trending_up = returns(trending_up);
returns_trending_down = returns(trending_down);

%% Moments

Moments_up(1,:) = mean(returns_trending_up);
Moments_up(2,:) = median(returns_trending_up);
Moments_up(3,:) = max(returns_trending_up);
Moments_up(4,:) = min(returns_trending_up);
Moments_up(5,:) = skewness(returns_trending_up);
Moments_up(6,:) = kurtosis(returns_trending_up);


Moments_down(1,:) = mean(returns_trending_down);
Moments_down(2,:) = median(returns_trending_down);
Moments_down(3,:) = max(returns_trending_down);
Moments_down(4,:) = min(returns_trending_down);
Moments_down(5,:) = skewness(returns_trending_down);
Moments_down(6,:) = kurtosis(returns_trending_down);

%% Generate normally distributed returns with variance of our sample - Just for comparison purpose

Sigma = nanstd(returns_all_one_column);
for i=1:size(returns_all_one_column,1)
    returns_normal(i,1)=normrnd(0,Sigma);
end

clearvars i;

%% Calculate amount of 2Sigma occurences



returns_trending_down((size(returns_trending_down,1)+1):(size(returns_all_one_column,1)))=NaN;
returns_trending_up((size(returns_trending_up,1)+1):(size(returns_all_one_column,1)))=NaN;
returns_comparison = [returns_all_one_column returns_normal returns_trending_up returns_trending_down];

Lower_sigma = returns_comparison < -Sigma;
Higher_sigma = returns_comparison > Sigma;

Amount_higher_sigma = nansum(Higher_sigma);
Amount_lower_sigma = nansum(Lower_sigma);

%% Calculate amount of 2Sigma occurences



returns_trending_down((size(returns_trending_down,1)+1):(size(returns_all_one_column,1)))=NaN;
returns_trending_up((size(returns_trending_up,1)+1):(size(returns_all_one_column,1)))=NaN;
returns_comparison = [returns_all_one_column returns_normal returns_trending_up returns_trending_down];

Lower_2sigma = returns_comparison < -2*Sigma;
Higher_2sigma = returns_comparison > 2*Sigma;

Amount_higher_2sigma = nansum(Higher_2sigma);
Amount_lower_2sigma = nansum(Lower_2sigma);



%% Calculate amount of 3Sigma occurences


Lower_3sigma = returns_comparison < -3*Sigma;
Higher_3sigma = returns_comparison > 3*Sigma;

Amount_higher_3sigma = nansum(Higher_3sigma);
Amount_lower_3sigma = nansum(Lower_3sigma);

clearvars Higher_2sigma
clearvars Higher_3sigma
clearvars Lower_2sigma
clearvars Lower_3sigma


%% Combine trending returns (We go short when trending down obvsiously)

returns_trending = vertcat(returns_trending_up, (-1*returns_trending_down));

%% Remove NaN values from the array
returns_trending(isnan(returns_trending))=[];



%% Bootstrap an equity curve based on the returns_trending
%% First we decide on how we will limit our downside. Every day we place a new SL to limit our downside as a multiple of the daily return std deviation
%% Then we randomly draw 2500 observations from our full sample of trending observations and we construct equity curves 100.000 times

limits = [-0.005 -0.01 -0.015]
for m=1:length(limits)
    
n=2500;
repl=100000;
returns_trending_random = returns_trending(randperm(length(returns_trending)));

%% Limit downside
limit = limits(m)
returns_trending_random(returns_trending_random<limit)=limit;
curves = zeros(n,repl);
curves(1,:) = 5000;

for k=1:repl
    index = round(rand(n,1)*n+0.5);
    returns_trending_bootstrap = zeros(n,1);
    for j=1:n
        returns_trending_bootstrap(j,1) = returns_trending_random(index(j,1),1);
    end
        for i=2:n
        curves(i,k) = curves(i-1,k).*(1+returns_trending_bootstrap(i-1,1));
        end
end

%% Then we construct some statistics per chosen limit level

final_equity = curves(n,:);
final_equity_avg(m) = mean(curves(n,:));
maxDD = maxdrawdown(curves);
maxDD_avg(m) = mean(maxDD);
maxDD_max(m) = max(maxDD);
end


clearvars i;
clearvars j;
clearvars k;
clearvars repl;
clearvars m