%% This progam will (hopefully) be a pairs trading strat
%% © Domenic Strazzulla 2009
%% Inputs

zscore = 60;
%this is the length of time to get the standard dev for z scoring

up =.9:.5:5;
down=.9:.5:5;
% up and down are how sensetive the trade signal is, the unit is Z
% scores. So a up or 3 and a down of 2 means that the signal will be hit if
% the differences in Z scores is more thatn 3 or less than - 2. (no need
% for a "-" sign on the down, it's taken care of later.) 


% The datasets are data and data2 - but you can change this to anything. 

%% Make Returns sequences

dataret1 = tick2ret(data);
% this is for 1 time unit

for i = 1:(size(data) - 5)
    dataret5(i,1) = ((data(i+5)-data(i))/data(i));
end
% this is for 5 time units

for i = 1:(size(data) - 10)
    dataret10(i,1) = ((data(i+10)-data(i))/data(i));
end
%10 units

for i = 1:(size(data) - 20)
    dataret20(i,1) = ((data(i+20)-data(i))/data(i));
end
%20 units


% Now we need standardize the 5, 10, and 20
for i = 1:(size(data) - zscore- 5) 
    zdata5(i,1) = dataret5(i+zscore)/std(dataret5(i:(zscore-1)+i));
end

for i = 1:(size(data) - zscore- 10) 
    zdata10(i,1) = dataret10(i+zscore)/std(dataret10(i:(zscore-1)+i));
end

for i = 1:(size(data) - zscore- 20) 
    zdata20(i,1) = dataret20(i+zscore)/std(dataret20(i:(zscore-1)+i));
end


% do it all again for the second dataset

data2ret1 = tick2ret(data2);
% this is for 1 time unit

for i = 1:(size(data2) - 5)
    data2ret5(i,1) = ((data2(i+5)-data2(i))/data2(i));
end
% this is for 5 time units

for i = 1:(size(data2) - 10)
    data2ret10(i,1) = ((data2(i+10)-data2(i))/data2(i));
end
%10 units

for i = 1:(size(data2) - 20)
    data2ret20(i,1) = ((data2(i+20)-data2(i))/data2(i));
end
%20 units

% Now we need standardize the 5, 10, and 20
for i = 1:(size(data2) - zscore- 5) 
    z2data5(i,1) = data2ret5(i+zscore)/std(data2ret5(i:(zscore-1)+i));
end

for i = 1:(size(data2) - zscore- 10) 
    z2data10(i,1) = data2ret10(i+zscore)/std(data2ret10(i:(zscore-1)+i));
end

for i = 1:(size(data2) - zscore- 20) 
    z2data20(i,1) = data2ret20(i+zscore)/std(data2ret20(i:(zscore-1)+i));
end

%% get the z-diffs
zdiff5 = (zdata5-z2data5);

zdiff10 = (zdata10-z2data10);

zdiff20 = (zdata20-z2data20);


%% Do a dicky fuller test 
dickyfullerdata = (dataret1- data2ret1);

[H,pValue,TestStat,CriticalValue] = dfARTest(dickyfullerdata);


%%
%plot the z scores
figure
plot(zdiff5);
xlabel('Number of Data Points')
ylabel('Z Score Diff 5')
title('5')

figure
plot(zdiff10);
xlabel('Number of Data Points')
ylabel('Z Score Diff 10')
title('10')

figure
plot(zdiff20);
xlabel('Number of Data Points')
ylabel('Z Score Diff 20')
title('20')

%% Now we have the basic framework to work from let's backtest to see if
%% this actually works


% This first one is for a z-diff 5 based strat, it's simple to change that
% to 10 or 20 though just replace 5 with 10 or 20

o=1;
t=1;

pos= zeros(size(zdiff5));

for d = up
    for s = down
        for i = 1 : size(zdiff5)
            if zdiff5(i)> d
                while zdiff5(i)  > 0
                    pos(i,1) = 1;
                    i= i +1;
                end
            elseif zdiff5(i) < -(s)
                while zdiff5(i) < 0
                    pos(i,1) = -1;
                    i = i +1;
                end
                %else
                %pos(i,1) = 0; don't use an else it messes things up(only 1
                %day long)
            end
        end
        for i = 1 : size(zdiff5)
            if pos(i) == 1
                ret(i,1) = (1 -(dataret1(i+zscore+5))) * (1 + (data2ret1(i+zscore+5)))- 1;
            elseif pos(i) == -1
                ret(i,1) = (1 +(dataret1(i+zscore+5))) * (1 - (data2ret1(i+zscore+5)))-1;
            else
                ret(i,1) = 0;
            end
        end
        %tempsize = size(zdiff5);
        sh = (sqrt(252)*mean(ret))/std(ret); % this becasue we are using daily data, if the interval was different we would have to change the 252.
        SH(o,t)=sh;
        pos= zeros(size(zdiff5));
        clear ret
        t=t+1;
    end
    o=o+1;
    t=1;
end

%The 3-d Sharpe Ratio Graph

SH2=SH; SH2(SH2==0)=NaN;
figure
surf(SH2), shading interp, lighting phong, 
xlabel('Down')
ylabel('Up')
zlabel('Sharpe')
view([80 35]), light('pos',[0.5, -0.9, 0.05])


%% © Domenic Strazzulla 2009