%% Important notes: 
% To run this script you need the dataset "Spot+DecFutures.xlsx".
% The code must be run sectionwise, as structured here. 
% Some of the subsections must be runned separately, such as plots and
% figures.
% If you have any issues with running the script, please contact
% sander.lie@student.bi.no

%% Part 1 - Preliminary analysis:
% Risk Premium statistics
clear;
clc;

% Import data and create variables:
data = readtable('Spot+DecFutures.xlsx');
data=data(~any(ismissing(data.Spot),2),:);
data(584:585,:) = [];
data = table2timetable(data);
data.RP19 = data.Fut19 - data.Spot;
data.RP18 = data.Fut18 - data.Spot;
data.RP17 = data.Fut17 - data.Spot;
data.RP16 = data.Fut16 - data.Spot;
data.RP15 = data.Fut15 - data.Spot;
data.RP14 = data.Fut14 - data.Spot;
data(1700:end,:) = [];

% Statistics:
Results = table();
Results.Observations(1) = 1699;
Results.Observations(2) = 1454;
Results.Observations(3) = 1207;
Results.Observations(4) = 983;
Results.Observations(5) = 729;
Results.Observations(6) = 485;

Results.Mean(1) = mean(data.RP19(1:1699));
Results.Mean(2) = mean(data.RP18(1:1454));
Results.Mean(3) = mean(data.RP17(1:1207));
Results.Mean(4) = mean(data.RP16(1:983));
Results.Mean(5) = mean(data.RP15(1:729));
Results.Mean(6) = mean(data.RP14(1:485));

Results.Std(1) = std(data.RP19(1:1699));
Results.Std(2) = std(data.RP18(1:1454));
Results.Std(3) = std(data.RP17(1:1207));
Results.Std(4) = std(data.RP16(1:983));
Results.Std(5) = std(data.RP15(1:729));
Results.Std(6) = std(data.RP14(1:485));

Results.Skew(1) = skewness(data.RP19(1:1699));
Results.Skew(2) = skewness(data.RP18(1:1454));
Results.Skew(3) = skewness(data.RP17(1:1207));
Results.Skew(4) = skewness(data.RP16(1:983));
Results.Skew(5) = skewness(data.RP15(1:729));
Results.Skew(6) = skewness(data.RP14(1:485));

Results.ExKurt(1) = kurtosis(data.RP19(1:1699))-3;
Results.ExKurt(2) = kurtosis(data.RP18(1:1454))-3;
Results.ExKurt(3) = kurtosis(data.RP17(1:1207))-3;
Results.ExKurt(4) = kurtosis(data.RP16(1:983))-3;
Results.ExKurt(5) = kurtosis(data.RP15(1:729))-3;
Results.ExKurt(6) = kurtosis(data.RP14(1:485))-3;

figure;
plot(data.Date, data.RP19)
hold on
plot(data.Date, data.RP18)
plot(data.Date, data.RP17)
plot(data.Date, data.RP16)
plot(data.Date, data.RP15)
plot(data.Date, data.RP14)
hold off
legend('RP 19', 'RP 18', 'RP 17', 'RP 16', 'RP 15', 'RP 14')
box('off')
xlabel('Year');
ylabel('Risk Premium in EUR')

% Historical prices of emission allowances
figure;
plot(data.Date,data.Spot)
hold on
plot(data.Date,data.Fut19);
plot(data.Date,data.Fut18);
plot(data.Date,data.Fut17);
plot(data.Date,data.Fut16);
plot(data.Date,data.Fut15);
plot(data.Date,data.Fut14);
legend('Spot', 'Fut19', 'Fut18', 'Fut17', 'Fut16', 'Fut15', 'Fut14', 'Location', 'SouthEast')
xlabel('Year');
ylabel('Price (in EUR)');
hold off

%% Part 2 - ADF tests
clear;
clc;

% Import data
data = readtable('Spot+DecFutures.xlsx');
data=data(~any(ismissing(data.Spot),2),:);
data(584:585,:) = [];
data = table2timetable(data);

%Creating logs:
data.Spot = log(data.Spot);
data.Fut19 = log(data.Fut19);
data.Fut18 = log(data.Fut18);
data.Fut17 = log(data.Fut17);
data.Fut16 = log(data.Fut16);
data.Fut15 = log(data.Fut15);
data.Fut14 = log(data.Fut14);
data.Electricity = log(data.Electricity);
data.Elec = data.Electricity;
data.Coal = log(data.Coal);
data.DAX = log(data.DAX);
data.NaturalGAS = log(data.NaturalGAS);
data.Ngas = data.NaturalGAS;

% Create a time series containing first differences
data.Spot_firstdiff = data.Spot - lagmatrix(data.Spot,1);
data.Fut19_firstdiff = data.Fut19 - lagmatrix(data.Fut19,1);
data.Fut18_firstdiff = data.Fut18 - lagmatrix(data.Fut18,1);
data.Fut17_firstdiff = data.Fut17 - lagmatrix(data.Fut17,1);
data.Fut16_firstdiff = data.Fut16 - lagmatrix(data.Fut16,1);
data.Fut15_firstdiff = data.Fut15 - lagmatrix(data.Fut15,1);
data.Fut14_firstdiff = data.Fut14 - lagmatrix(data.Fut14,1);
data.Oil_fd = data.Oil - lagmatrix(data.Oil,1);
data.Elec_fd = data.Elec - lagmatrix(data.Elec,1);
data.Coal_fd = data.Coal - lagmatrix(data.Coal,1);
data.DAX_fd = data.DAX - lagmatrix(data.DAX,1);
data.Ngas_fd = data.Ngas - lagmatrix(data.Ngas,1);

% ADF Test (Unit root test) - Check each time series for stationarity
results = table();

i = 0:2;
for model = i
    maxlag = 12;
    ic = 'AIC';
    alpha = [0.01; 0.05; 0.10];
    model
    

    % Spot:

    % Levels
    [adfstat, pval, critval, ~, lags_level] = augdfautolag(data.Spot, model, maxlag, ic);
    disp('Spot Levels');
    disp(table(pval, adfstat, lags_level));
    critval = critval(1:3);
    reject = adfstat > critval;
    disp(table(alpha, critval, reject));
    if model == 0
        results.SpotLevels(1) = adfstat;
    end
    if model == 1
        results.SpotLevels(2) = adfstat;
    end
    if model == 2
        results.SpotLevels(3) = adfstat;
    end


    % First differences
    [adfstat, pval, critval, ~, lags_firstdiff] = augdfautolag(data.Spot_firstdiff(2:end), model, maxlag, ic);
    disp('Spot first differences');
    disp(table(pval, adfstat, lags_firstdiff));
    critval = critval(1:3);
    reject = adfstat > critval;
    disp(table(alpha, critval, reject));
    if model == 0
        results.SpotDiff(1) = adfstat;
    end
    if model == 1
        results.SpotDiff(2) = adfstat;
    end
    if model == 2
        results.SpotDiff(3) = adfstat;
    end
    
    % Dec2019 Future
    % Levels
    [adfstat, pval, critval, ~, lags_level] = augdfautolag(data.Fut19(1:1699), model, maxlag, ic);
    disp('Fut19 Levels');
    disp(table(pval, adfstat, lags_level));
    critval = critval(1:3);
    reject = adfstat > critval;
    disp(table(alpha, critval, reject));
    if model == 0
        results.Fut19Levels(1) = adfstat;
    end
    if model == 1
        results.Fut19Levels(2) = adfstat;
    end
    if model == 2
        results.Fut19Levels(3) = adfstat;
    end

    % First differences
    [adfstat, pval, critval, ~, lags_firstdiff] = augdfautolag(data.Fut19_firstdiff(2:1699), model, maxlag, ic);
    disp('Fut19 first differences');
    disp(table(pval, adfstat, lags_firstdiff));
    critval = critval(1:3);
    reject = adfstat > critval;
    disp(table(alpha, critval, reject));
    if model == 0
        results.Fut19Diff(1) = adfstat;
    end
    if model == 1
        results.Fut19Diff(2) = adfstat;
    end
    if model == 2
        results.Fut19Diff(3) = adfstat;
    end
    
    
    % Dec2018 Future
    % Levels
    [adfstat, pval, critval, ~, lags_level] = augdfautolag(data.Fut18(1:1454), model, maxlag, ic);
    disp('Fut18 Levels');
    disp(table(pval, adfstat, lags_level));
    critval = critval(1:3);
    reject = adfstat > critval;
    disp(table(alpha, critval, reject));
    if model == 0
        results.Fut18Levels(1) = adfstat;
    end
    if model == 1
        results.Fut18Levels(2) = adfstat;
    end
    if model == 2
        results.Fut18Levels(3) = adfstat;
    end
    
    % First differences
    [adfstat, pval, critval, ~, lags_firstdiff] = augdfautolag(data.Fut18_firstdiff(2:1454), model, maxlag, ic);
    disp('Fut18 first differences');
    disp(table(pval, adfstat, lags_firstdiff));
    critval = critval(1:3);
    reject = adfstat > critval;
    disp(table(alpha, critval, reject));
    if model == 0
        results.Fut18Diff(1) = adfstat;
    end
    if model == 1
        results.Fut18Diff(2) = adfstat;
    end
    if model == 2
        results.Fut18Diff(3) = adfstat;
    end
    

    % Dec2017 Future
    % Levels
    [adfstat, pval, critval, ~, lags_level] = augdfautolag(data.Fut17(1:1207), model, maxlag, ic);
    disp('Fut17 Levels');
    disp(table(pval, adfstat, lags_level));
    critval = critval(1:3);
    reject = adfstat > critval;
    disp(table(alpha, critval, reject));
    if model == 0
        results.Fut17Levels(1) = adfstat;
    end
    if model == 1
        results.Fut17Levels(2) = adfstat;
    end
    if model == 2
        results.Fut17Levels(3) = adfstat;
    end
    
    % First differences
    [adfstat, pval, critval, ~, lags_firstdiff] = augdfautolag(data.Fut17_firstdiff(2:1207), model, maxlag, ic);
    disp('Fut17 first differences');
    disp(table(pval, adfstat, lags_firstdiff));
    critval = critval(1:3);
    reject = adfstat > critval;
    disp(table(alpha, critval, reject));
    if model == 0
        results.Fut17Diff(1) = adfstat;
    end
    if model == 1
        results.Fut17Diff(2) = adfstat;
    end
    if model == 2
        results.Fut17Diff(3) = adfstat;
    end
    
    
    % Dec2016 Future
    % Levels
    [adfstat, pval, critval, ~, lags_level] = augdfautolag(data.Fut16(1:983), model, maxlag, ic);
    disp('Fut16 Levels');
    disp(table(pval, adfstat, lags_level));
    critval = critval(1:3);
    reject = adfstat > critval;
    disp(table(alpha, critval, reject));
    if model == 0
        results.Fut16Levels(1) = adfstat;
    end
    if model == 1
        results.Fut16Levels(2) = adfstat;
    end
    if model == 2
        results.Fut16Levels(3) = adfstat;
    end
    
    % First differences
    [adfstat, pval, critval, ~, lags_firstdiff] = augdfautolag(data.Fut16_firstdiff(2:983), model, maxlag, ic);
    disp('Fut16 first differences');
    disp(table(pval, adfstat, lags_firstdiff));
    critval = critval(1:3);
    reject = adfstat > critval;
    disp(table(alpha, critval, reject));
    if model == 0
        results.Fut16Diff(1) = adfstat;
    end
    if model == 1
        results.Fut16Diff(2) = adfstat;
    end
    if model == 2
        results.Fut16Diff(3) = adfstat;
    end
    
    % Dec2015 Future
    % Levels
    [adfstat, pval, critval, ~, lags_level] = augdfautolag(data.Fut15(1:729), model, maxlag, ic);
    disp('Fut15 Levels');
    disp(table(pval, adfstat, lags_level));
    critval = critval(1:3);
    reject = adfstat > critval;
    disp(table(alpha, critval, reject));
    if model == 0
        results.Fut15Levels(1) = adfstat;
    end
    if model == 1
        results.Fut15Levels(2) = adfstat;
    end
    if model == 2
        results.Fut15Levels(3) = adfstat;
    end

    % First differences
    [adfstat, pval, critval, ~, lags_firstdiff] = augdfautolag(data.Fut15_firstdiff(2:729), model, maxlag, ic);
    disp('Fut15 first differences');
    disp(table(pval, adfstat, lags_firstdiff));
    critval = critval(1:3);
    reject = adfstat > critval;
    disp(table(alpha, critval, reject));
    if model == 0
        results.Fut15Diff(1) = adfstat;
    end
    if model == 1
        results.Fut15Diff(2) = adfstat;
    end
    if model == 2
        results.Fut15Diff(3) = adfstat;
    end

    % Dec2014 Future
    % Levels
    [adfstat, pval, critval, ~, lags_level] = augdfautolag(data.Fut14 (1:485), model, maxlag, ic);
    disp('Fut14 Levels');
    disp(table(pval, adfstat, lags_level));
    critval = critval(1:3);
    reject = adfstat > critval;
    disp(table(alpha, critval, reject));
    if model == 0
        results.Fut14Levels(1) = adfstat;
    end
    if model == 1
        results.Fut14Levels(2) = adfstat;
    end
    if model == 2
        results.Fut14Levels(3) = adfstat;
    end

    % First differences
    [adfstat, pval, critval, ~, lags_firstdiff] = augdfautolag(data.Fut14_firstdiff(2:485), model, maxlag, ic);
    disp('Fut14 first differences');
    disp(table(pval, adfstat, lags_firstdiff));
    critval = critval(1:3);
    reject = adfstat > critval;
    disp(table(alpha, critval, reject));
    if model == 0
        results.Fut14Diff(1) = adfstat;
    end
    if model == 1
        results.Fut14Diff(2) = adfstat;
    end
    if model == 2
        results.Fut14Diff(3) = adfstat;
    end
    
    % Oil:
    % Levels
    [adfstat, pval, critval, ~, lags_level] = augdfautolag(data.Oil, model, maxlag, ic);
    disp('Oil Levels');
    disp(table(pval, adfstat, lags_level));
    critval = critval(1:3);
    reject = adfstat > critval;
    disp(table(alpha, critval, reject));
    if model == 0
        results.OilLevels(1) = adfstat;
    end
    if model == 1
        results.OilLevels(2) = adfstat;
    end
    if model == 2
        results.OilLevels(3) = adfstat;
    end


    % First differences
    [adfstat, pval, critval, ~, lags_firstdiff] = augdfautolag(data.Oil_fd(2:end), model, maxlag, ic);
    disp('Oil first differences');
    disp(table(pval, adfstat, lags_firstdiff));
    critval = critval(1:3);
    reject = adfstat > critval;
    disp(table(alpha, critval, reject));
    if model == 0
        results.OilDiff(1) = adfstat;
    end
    if model == 1
        results.OilDiff(2) = adfstat;
    end
    if model == 2
        results.OilDiff(3) = adfstat;
    end
    
    % Elec:
    % Levels
    [adfstat, pval, critval, ~, lags_level] = augdfautolag(data.Elec, model, maxlag, ic);
    disp('Elec Levels');
    disp(table(pval, adfstat, lags_level));
    critval = critval(1:3);
    reject = adfstat > critval;
    disp(table(alpha, critval, reject));
    if model == 0
        results.ElecLevels(1) = adfstat;
    end
    if model == 1
        results.ElecLevels(2) = adfstat;
    end
    if model == 2
        results.ElecLevels(3) = adfstat;
    end

    % First differences
    [adfstat, pval, critval, ~, lags_firstdiff] = augdfautolag(data.Elec_fd(2:end), model, maxlag, ic);
    disp('Elec first differences');
    disp(table(pval, adfstat, lags_firstdiff));
    critval = critval(1:3);
    reject = adfstat > critval;
    disp(table(alpha, critval, reject));
    if model == 0
        results.ElecDiff(1) = adfstat;
    end
    if model == 1
        results.ElecDiff(2) = adfstat;
    end
    if model == 2
        results.ElecDiff(3) = adfstat;
    end
    
    
    % Coal:
    % Levels
    [adfstat, pval, critval, ~, lags_level] = augdfautolag(data.Coal(557:end), model, maxlag, ic);
    disp('Coal Levels');
    disp(table(pval, adfstat, lags_level));
    critval = critval(1:3);
    reject = adfstat > critval;
    disp(table(alpha, critval, reject));
    if model == 0
        results.CoalLevels(1) = adfstat;
    end
    if model == 1
        results.CoalLevels(2) = adfstat;
    end
    if model == 2
        results.CoalLevels(3) = adfstat;
    end

    % First differences
    [adfstat, pval, critval, ~, lags_firstdiff] = augdfautolag(data.Coal_fd(558:end), model, maxlag, ic);
    disp('Coal first differences');
    disp(table(pval, adfstat, lags_firstdiff));
    critval = critval(1:3);
    reject = adfstat > critval;
    disp(table(alpha, critval, reject));
    if model == 0
        results.CoalDiff(1) = adfstat;
    end
    if model == 1
        results.CoalDiff(2) = adfstat;
    end
    if model == 2
        results.CoalDiff(3) = adfstat;
    end
    
    % DAX:
    % DAX data:
    DAXdata = table();
    DAXdata.DAX = data.DAX;
    DAXdata.DAX_fd = data.DAX_fd;
    DAXdata=DAXdata(~any(ismissing(DAXdata.DAX),2),:);
    DAXdata=DAXdata(~any(ismissing(DAXdata.DAX_fd),2),:);
    
    % Levels
    [adfstat, pval, critval, ~, lags_level] = augdfautolag(DAXdata.DAX(1:end), model, maxlag, ic);
    disp('DAX Levels');
    disp(table(pval, adfstat, lags_level));
    critval = critval(1:3);
    reject = adfstat > critval;
    disp(table(alpha, critval, reject));
    if model == 0
        results.DAXLevels(1) = adfstat;
    end
    if model == 1
        results.DAXLevels(2) = adfstat;
    end
    if model == 2
        results.DAXLevels(3) = adfstat;
    end

    % First differences
    [adfstat, pval, critval, ~, lags_firstdiff] = augdfautolag(DAXdata.DAX_fd(1:end), model, maxlag, ic);
    disp('DAX first differences');
    disp(table(pval, adfstat, lags_firstdiff));
    critval = critval(1:3);
    reject = adfstat > critval;
    disp(table(alpha, critval, reject));
    if model == 0
        results.DAXDiff(1) = adfstat;
    end
    if model == 1
        results.DAXDiff(2) = adfstat;
    end
    if model == 2
        results.DAXDiff(3) = adfstat;
    end
    
    % Ngas:
    % Ngas data:
    Ngasdata = table();
    Ngasdata.Ngas = data.Ngas;
    Ngasdata.Ngas_fd = data.Ngas_fd;
    Ngasdata=Ngasdata(~any(ismissing(Ngasdata.Ngas),2),:);
    Ngasdata=Ngasdata(~any(ismissing(Ngasdata.Ngas_fd),2),:);
    
    % Levels
    [adfstat, pval, critval, ~, lags_level] = augdfautolag(Ngasdata.Ngas(1:end), model, maxlag, ic);
    disp('Ngas Levels');
    disp(table(pval, adfstat, lags_level));
    critval = critval(1:3);
    reject = adfstat > critval;
    disp(table(alpha, critval, reject));
    if model == 0
        results.NgasLevels(1) = adfstat;
    end
    if model == 1
        results.NgasLevels(2) = adfstat;
    end
    if model == 2
        results.NgasLevels(3) = adfstat;
    end

    % First differences
    [adfstat, pval, critval, ~, lags_firstdiff] = augdfautolag(Ngasdata.Ngas_fd(1:end), model, maxlag, ic);
    disp('Ngas first differences');
    disp(table(pval, adfstat, lags_firstdiff));
    critval = critval(1:3);
    reject = adfstat > critval;
    disp(table(alpha, critval, reject));
    if model == 0
        results.NgasDiff(1) = adfstat;
    end
    if model == 1
        results.NgasDiff(2) = adfstat;
    end
    if model == 2
        results.NgasDiff(3) = adfstat;
    end
end
%% Part 3 - Cointegration (Johansen test):
clear;
clc;

% Import data and create variables:
data = readtable('Spot+DecFutures.xlsx');
data=data(~any(ismissing(data.Spot),2),:);
data(584:585,:) = [];
data = table2timetable(data);
data.Spot = log(data.Spot);
data.Fut19 = log(data.Fut19);
data.Fut18 = log(data.Fut18);
data.Fut17 = log(data.Fut17);
data.Fut16 = log(data.Fut16);
data.Fut15 = log(data.Fut15);
data.Fut14 = log(data.Fut14);

% Determining lag levels:
% Spot and Fut19:
reg1 = table();
reg1.Spot = data.Spot;
reg1.Fut19 = data.Fut19;
reg1(1700:end,:) = [];
reg1 = reg1{:,:};
[laglength, AIC, logL] = VARlag(reg1,12)
%Optimal = 12

% Spot and Fut18:
reg2 = table();
reg2.Spot = data.Spot;
reg2.Fut18 = data.Fut18;
reg2(1455:end,:) = [];
reg2 = reg2{:,:};
[laglength, AIC, logL] = VARlag(reg2,12)
%Optimal = 8

% Spot and Fut17:
reg3 = table();
reg3.Spot = data.Spot;
reg3.Fut17 = data.Fut17;
reg3(1208:end,:) = [];
reg3 = reg3{:,:};
[laglength, AIC, logL] = VARlag(reg3,12)
%Optimal = 12

%Spot and Fut16
reg4 = table();
reg4.Spot = data.Spot;
reg4.Fut16 = data.Fut16;
reg4(984:end,:) = [];
reg4 = reg4{:,:};
[laglength, AIC, logL] = VARlag(reg4,12)
%Optimal = 7

% Spot and Fut15
reg5 = table();
reg5.Spot = data.Spot;
reg5.Fut15 = data.Fut15;
reg5(730:end,:) = [];
reg5 = reg5{:,:};
[laglength, AIC, logL] = VARlag(reg5,12)
%Optimal = 10

%Spot and Fut14
reg6 = table();
reg6.Spot = data.Spot;
reg6.Fut14 = data.Fut14;
reg6(486:end,:) = [];
reg6 = reg6{:,:};
[laglength, AIC, logL] = VARlag(reg6,12)
%Optimal = 7

%Tables for results:
traceH0 = table();
MaxEigH0 = table();
traceH1 = table();
MaxEigH1 = table();

% Spot and Fut19:
% Null hypothesis - cointegrated series with intercept but no trend
[~,pValue,~,~,~] = jcitest(reg1, 'model', 'H1*','lags', 12, 'test', 'trace', 'display', 'off'); 
disp('trace H0'); 
disp(pValue);
traceH0.Fut19(1) = pValue.r0;
traceH0.Fut19(2) = pValue.r1;

[~,pValue,~,~,~] = jcitest(reg1, 'model', 'H1*','lags', 12, 'test', 'maxeig', 'display', 'off'); 
disp('max eigenvector H0'); 
disp(pValue);
MaxEigH0.Fut19(1) = pValue.r0;
MaxEigH0.Fut19(2) = pValue.r1;

% Null hypothesis - cointegrated series with no intercept, no trend
[~,pValue,~,~,~] = jcitest(reg1, 'model', 'H2','lags', 12, 'test', 'trace', 'display', 'off'); 
disp('trace H1');
disp(pValue)
traceH1.Fut19(1) = pValue.r0;
traceH1.Fut19(2) = pValue.r1;

[~,pValue,~,~,~] = jcitest(reg1, 'model', 'H2','lags', 12, 'test', 'maxeig', 'display', 'off'); 
disp('max eigenvector H1'); 
disp(pValue);
MaxEigH1.Fut19(1) = pValue.r0;
MaxEigH1.Fut19(2) = pValue.r1;

% Spot and Fut18:
% Null hypothesis - cointegrated series with intercept but no trend
[~,pValue,~,~,~] = jcitest(reg2, 'model', 'H1*','lags', 8, 'test', 'trace', 'display', 'off'); 
disp('trace H0'); 
disp(pValue);
traceH0.Fut18(1) = pValue.r0;
traceH0.Fut18(2) = pValue.r1;

[~,pValue,~,~,~] = jcitest(reg2, 'model', 'H1*','lags', 8, 'test', 'maxeig', 'display', 'off'); 
disp('max eigenvector H0'); 
disp(pValue);
MaxEigH0.Fut18(1) = pValue.r0;
MaxEigH0.Fut18(2) = pValue.r1;

% Null hypothesis - cointegrated series with no intercept, no trend
[~,pValue,~,~,~] = jcitest(reg2, 'model', 'H2','lags', 8, 'test', 'trace', 'display', 'off'); 
disp('trace H1');
disp(pValue);
traceH1.Fut18(1) = pValue.r0;
traceH1.Fut18(2) = pValue.r1;

[~,pValue,~,~,~] = jcitest(reg2, 'model', 'H2','lags', 8, 'test', 'maxeig', 'display', 'off'); 
disp('max eigenvector H1'); 
disp(pValue);
MaxEigH1.Fut18(1) = pValue.r0;
MaxEigH1.Fut18(2) = pValue.r1;

% Spot and Fut17:
% Null hypothesis - cointegrated series with intercept but no trend
[~,pValue,~,~,~] = jcitest(reg3, 'model', 'H1*','lags', 12, 'test', 'trace', 'display', 'off'); 
disp('trace H0'); 
disp(pValue);
traceH0.Fut17(1) = pValue.r0;
traceH0.Fut17(2) = pValue.r1;

[~,pValue,~,~,~] = jcitest(reg3, 'model', 'H1*','lags', 12, 'test', 'maxeig', 'display', 'off'); 
disp('max eigenvector H0'); 
disp(pValue);
MaxEigH0.Fut17(1) = pValue.r0;
MaxEigH0.Fut17(2) = pValue.r1;

% Null hypothesis - cointegrated series with no intercept, no trend
[~,pValue,~,~,~] = jcitest(reg3, 'model', 'H2','lags', 12, 'test', 'trace', 'display', 'off'); 
disp('trace H1');
disp(pValue);
traceH1.Fut17(1) = pValue.r0;
traceH1.Fut17(2) = pValue.r1;

[~,pValue,~,~,~] = jcitest(reg3, 'model', 'H2','lags', 12, 'test', 'maxeig', 'display', 'off'); 
disp('max eigenvector H1'); 
disp(pValue);
MaxEigH1.Fut17(1) = pValue.r0;
MaxEigH1.Fut17(2) = pValue.r1;

% Spot and Fut16:
% Null hypothesis - cointegrated series with intercept but no trend
[~,pValue,~,~,~] = jcitest(reg4, 'model', 'H1*','lags', 7, 'test', 'trace', 'display', 'off'); 
disp('trace H0'); 
disp(pValue);
traceH0.Fut16(1) = pValue.r0;
traceH0.Fut16(2) = pValue.r1;

[~,pValue,~,~,~] = jcitest(reg4, 'model', 'H1*','lags', 7, 'test', 'maxeig', 'display', 'off'); 
disp('max eigenvector H0'); 
disp(pValue);
MaxEigH0.Fut16(1) = pValue.r0;
MaxEigH0.Fut16(2) = pValue.r1;

% Null hypothesis - cointegrated series with no intercept, no trend
[~,pValue,~,~,~] = jcitest(reg4, 'model', 'H2','lags', 7, 'test', 'trace', 'display', 'off'); 
disp('trace H1');
disp(pValue);
traceH1.Fut16(1) = pValue.r0;
traceH1.Fut16(2) = pValue.r1;

[~,pValue,~,~,~] = jcitest(reg4, 'model', 'H2','lags', 7, 'test', 'maxeig', 'display', 'off'); 
disp('max eigenvector H1'); 
disp(pValue);
MaxEigH1.Fut16(1) = pValue.r0;
MaxEigH1.Fut16(2) = pValue.r1;

% Spot and Fut15:
% Null hypothesis - cointegrated series with intercept but no trend
[~,pValue,~,~,~] = jcitest(reg5, 'model', 'H1*','lags', 10, 'test', 'trace', 'display', 'off'); 
disp('trace H0'); 
disp(pValue);
traceH0.Fut15(1) = pValue.r0;
traceH0.Fut15(2) = pValue.r1;

[~,pValue,~,~,~] = jcitest(reg5, 'model', 'H1*','lags', 10, 'test', 'maxeig', 'display', 'off'); 
disp('max eigenvector H0'); 
disp(pValue);
MaxEigH0.Fut15(1) = pValue.r0;
MaxEigH0.Fut15(2) = pValue.r1;

% Null hypothesis - cointegrated series with no intercept, no trend
[~,pValue,~,~,~] = jcitest(reg5, 'model', 'H2','lags', 10, 'test', 'trace', 'display', 'off'); 
disp('trace H1');
disp(pValue);
traceH1.Fut15(1) = pValue.r0;
traceH1.Fut15(2) = pValue.r1;

[~,pValue,~,~,~] = jcitest(reg5, 'model', 'H2','lags', 10, 'test', 'maxeig', 'display', 'off'); 
disp('max eigenvector H1'); 
disp(pValue);
MaxEigH1.Fut15(1) = pValue.r0;
MaxEigH1.Fut15(2) = pValue.r1;

% Spot and Fut14:
% Null hypothesis - cointegrated series with intercept but no trend
[~,pValue,~,~,~] = jcitest(reg6, 'model', 'H1*','lags', 7, 'test', 'trace', 'display', 'off'); 
disp('trace H0'); 
disp(pValue);
traceH0.Fut14(1) = pValue.r0;
traceH0.Fut14(2) = pValue.r1;

[~,pValue,~,~,~] = jcitest(reg6, 'model', 'H1*','lags', 7, 'test', 'maxeig', 'display', 'off'); 
disp('max eigenvector H0'); 
disp(pValue);
MaxEigH0.Fut14(1) = pValue.r0;
MaxEigH0.Fut14(2) = pValue.r1;

% Null hypothesis - cointegrated series with no intercept, no trend
[~,pValue,~,~,~] = jcitest(reg6, 'model', 'H2','lags', 7, 'test', 'trace', 'display', 'off'); 
disp('trace H1');
disp(pValue);
traceH1.Fut14(1) = pValue.r0;
traceH1.Fut14(2) = pValue.r1;

[~,pValue,~,~,~] = jcitest(reg6, 'model', 'H2','lags', 7, 'test', 'maxeig', 'display', 'off'); 
disp('max eigenvector H1'); 
disp(pValue);
MaxEigH1.Fut14(1) = pValue.r0;
MaxEigH1.Fut14(2) = pValue.r1;

%% Part 4 - ECM and ECM-GARCH:
clear;
clc;

% Import data and create variables:
data = readtable('Spot+DecFutures.xlsx');
data=data(~any(ismissing(data.Spot),2),:);
data(584:585,:) = [];
data = table2timetable(data);
data.Spot = log(data.Spot);
data.Fut19 = log(data.Fut19);
data.Fut18 = log(data.Fut18);
data.Fut17 = log(data.Fut17);
data.Fut16 = log(data.Fut16);
data.Fut15 = log(data.Fut15);
data.Fut14 = log(data.Fut14);

% Creating first differences:
data.Spot_firstdiff = data.Spot - lagmatrix(data.Spot,1);
data.Fut19_firstdiff = data.Fut19 - lagmatrix(data.Fut19,1);
data.Fut18_firstdiff = data.Fut18 - lagmatrix(data.Fut18,1);
data.Fut17_firstdiff = data.Fut17 - lagmatrix(data.Fut17,1);
data.Fut16_firstdiff = data.Fut16 - lagmatrix(data.Fut16,1);
data.Fut15_firstdiff = data.Fut15 - lagmatrix(data.Fut15,1);
data.Fut14_firstdiff = data.Fut14 - lagmatrix(data.Fut14,1);

% Finding residuals between Spot and lagged Fut19:
reg1 = table();
reg1.Date = data.Date;
reg1.Spot = data.Spot;
reg1.LagFut19 = lagmatrix(data.Fut19,1);
reg1(1,:) = [];
reg1(1700:end,:) = [];
model1 = fitlm(reg1, 'Spot ~ LagFut19')
reg1.residuals1 = model1.Residuals.Raw;
figure;
plot(reg1.Date, reg1.residuals1);

% Defining ECM variables - Spot and Fut19:
ECM1 = table();
ECM1.Date = reg1.Date;
ECM1.Spot_firstdiff = data.Spot_firstdiff(2:1700);
ECM1.LaggedResiduals = lagmatrix(reg1.residuals1,1);
%Lag variables - Fut:
ECM1.Fut19_firstdiff_lag1 = reg1.LagFut19 - lagmatrix(reg1.LagFut19,1);
ECM1.Fut19_firstdiff_lag2 = lagmatrix(ECM1.Fut19_firstdiff_lag1,1);
ECM1.Fut19_firstdiff_lag3 = lagmatrix(ECM1.Fut19_firstdiff_lag1,2);
ECM1.Fut19_firstdiff_lag4 = lagmatrix(ECM1.Fut19_firstdiff_lag1,3);
ECM1.Fut19_firstdiff_lag5 = lagmatrix(ECM1.Fut19_firstdiff_lag1,4);
ECM1.Fut19_firstdiff_lag6 = lagmatrix(ECM1.Fut19_firstdiff_lag1,5);
ECM1.Fut19_firstdiff_lag7 = lagmatrix(ECM1.Fut19_firstdiff_lag1,6);
ECM1.Fut19_firstdiff_lag8 = lagmatrix(ECM1.Fut19_firstdiff_lag1,7);
ECM1.Fut19_firstdiff_lag9 = lagmatrix(ECM1.Fut19_firstdiff_lag1,8);
ECM1.Fut19_firstdiff_lag10 = lagmatrix(ECM1.Fut19_firstdiff_lag1,9);
ECM1.Fut19_firstdiff_lag11 = lagmatrix(ECM1.Fut19_firstdiff_lag1,10);
ECM1.Fut19_firstdiff_lag12 = lagmatrix(ECM1.Fut19_firstdiff_lag1,11);
%Lag variables spot:
ECM1.Spot_firstdiff_lag1 = lagmatrix(ECM1.Spot_firstdiff,1);
ECM1.Spot_firstdiff_lag2 = lagmatrix(ECM1.Spot_firstdiff,2);
ECM1.Spot_firstdiff_lag3 = lagmatrix(ECM1.Spot_firstdiff,3);
ECM1.Spot_firstdiff_lag4 = lagmatrix(ECM1.Spot_firstdiff,4);
ECM1.Spot_firstdiff_lag5 = lagmatrix(ECM1.Spot_firstdiff,5);
ECM1.Spot_firstdiff_lag6 = lagmatrix(ECM1.Spot_firstdiff,6);
ECM1.Spot_firstdiff_lag7 = lagmatrix(ECM1.Spot_firstdiff,7);
ECM1.Spot_firstdiff_lag8 = lagmatrix(ECM1.Spot_firstdiff,8);
ECM1.Spot_firstdiff_lag9 = lagmatrix(ECM1.Spot_firstdiff,9);
ECM1.Spot_firstdiff_lag10 = lagmatrix(ECM1.Spot_firstdiff,10);
ECM1.Spot_firstdiff_lag11 = lagmatrix(ECM1.Spot_firstdiff,11);
ECM1.Spot_firstdiff_lag12 = lagmatrix(ECM1.Spot_firstdiff,12);

ECM_1 = fitlm(ECM1, ['Spot_firstdiff ~ LaggedResiduals + Fut19_firstdiff_lag1 + Fut19_firstdiff_lag2' ... 
    '+ Fut19_firstdiff_lag3 + Fut19_firstdiff_lag4 + Fut19_firstdiff_lag5 + Fut19_firstdiff_lag6' ...
    '+ Fut19_firstdiff_lag7 + Fut19_firstdiff_lag8 + Fut19_firstdiff_lag9 + Fut19_firstdiff_lag10' ...
    '+ Fut19_firstdiff_lag11 + Fut19_firstdiff_lag12 + Spot_firstdiff_lag1 + Spot_firstdiff_lag2' ...
    '+ Spot_firstdiff_lag3 + Spot_firstdiff_lag4 + Spot_firstdiff_lag5 + Spot_firstdiff_lag6' ...
    '+ Spot_firstdiff_lag7 + Spot_firstdiff_lag8 + Spot_firstdiff_lag9 + Spot_firstdiff_lag10' ...
    '+ Spot_firstdiff_lag11 + Spot_firstdiff_lag12'])

% Engle's ARCH Test for ECM_F19:
ECM1.SqResiduals = (ECM_1.Residuals.Raw).^2;
ECM1.res1 = lagmatrix(ECM1.SqResiduals,1);
ECM1.res2 = lagmatrix(ECM1.SqResiduals,2);
ECM1.res3 = lagmatrix(ECM1.SqResiduals,3);
ECM1.res4 = lagmatrix(ECM1.SqResiduals,4);
ECM1.res5 = lagmatrix(ECM1.SqResiduals,5);
ECM1.res6 = lagmatrix(ECM1.SqResiduals,6);
ECM1.res7 = lagmatrix(ECM1.SqResiduals,7);
ECM1.res8 = lagmatrix(ECM1.SqResiduals,8);
ECM1.res9 = lagmatrix(ECM1.SqResiduals,9);
ECM1.res10 = lagmatrix(ECM1.SqResiduals,10);
ECM1.res11 = lagmatrix(ECM1.SqResiduals,11);
ECM1.res12 = lagmatrix(ECM1.SqResiduals,12);

ECM1 = fillmissing(ECM1, 'constant', 0, 'DataVariables', ...
                    {'res1','res2','res3','res4','res5','res6','res7', ...
                    'res8','res9','res10','res11','res12'});

% Running the auxiliary regression
AuxReg1 = fitlm(ECM1, ['SqResiduals ~ res1 + res2 + res3 + res3 + res4 + res5' ...
    '+ res6 + res7 + res8 + res9 + res10 + res11 + res12'])

% Testing for autocorrelation:
T =  1687;
alpha = [0.1, 0.05, 0.01]
TestStatARCH = (T)*AuxReg1.Rsquared.Ordinary
pVal = 1 - cdf('Chisquare', TestStatARCH, 12)
criticalARCH = chi2inv(1-alpha, 12);
rejection = TestStatARCH > criticalARCH

% ECM-GARCH: Spot and Fut19:
GARCH1 = garch('GARCHlags', 1, 'ARCHlags', 1)
[estMdl, estParamCov, logL] = estimate(GARCH1, ECM1.Spot_firstdiff)
condVar = infer(estMdl, ECM1.Spot_firstdiff);
ECM1.condVol = sqrt(condVar);

ECM_GARCH_1 = fitlm(ECM1, ['Spot_firstdiff ~ LaggedResiduals + Fut19_firstdiff_lag1 + Fut19_firstdiff_lag2' ... 
    '+ Fut19_firstdiff_lag3 + Fut19_firstdiff_lag4 + Fut19_firstdiff_lag5 + Fut19_firstdiff_lag6' ...
    '+ Fut19_firstdiff_lag7 + Fut19_firstdiff_lag8 + Fut19_firstdiff_lag9 + Fut19_firstdiff_lag10' ...
    '+ Fut19_firstdiff_lag11 + Fut19_firstdiff_lag12 + Spot_firstdiff_lag1 + Spot_firstdiff_lag2' ...
    '+ Spot_firstdiff_lag3 + Spot_firstdiff_lag4 + Spot_firstdiff_lag5 + Spot_firstdiff_lag6' ...
    '+ Spot_firstdiff_lag7 + Spot_firstdiff_lag8 + Spot_firstdiff_lag9 + Spot_firstdiff_lag10' ...
    '+ Spot_firstdiff_lag11 + Spot_firstdiff_lag12 + condVol'])

figure;
plot(ECM1.Date, ECM1.Spot_firstdiff);
hold on;
plot(ECM1.Date, ECM1.condVol);
hold off;
legend('Log-differenced spot price', 'Inferred volatility');
box('off');

% Finding residuals between Spot and lagged Fut18:
reg2 = table();
reg2.Date = data.Date;
reg2.Spot = data.Spot;
reg2.LagFut18 = lagmatrix(data.Fut18,1);
reg2(1,:) = [];
reg2(1455:end,:) = [];
model2 = fitlm(reg2, 'Spot ~ LagFut18')
reg2.residuals = model2.Residuals.Raw;
plot(reg2.Date, reg2.residuals)

% Defining ECM variables - Spot and Fut18:
ECM2 = table();
ECM2.Date = reg2.Date;
ECM2.Spot_firstdiff = data.Spot_firstdiff(2:1455);
ECM2.LaggedResiduals = lagmatrix(reg2.residuals,1);
%Lag variables - Fut:
ECM2.Fut18_firstdiff_lag1 = reg2.LagFut18 - lagmatrix(reg2.LagFut18,1);
ECM2.Fut18_firstdiff_lag2 = lagmatrix(ECM2.Fut18_firstdiff_lag1,1);
ECM2.Fut18_firstdiff_lag3 = lagmatrix(ECM2.Fut18_firstdiff_lag1,2);
ECM2.Fut18_firstdiff_lag4 = lagmatrix(ECM2.Fut18_firstdiff_lag1,3);
ECM2.Fut18_firstdiff_lag5 = lagmatrix(ECM2.Fut18_firstdiff_lag1,4);
ECM2.Fut18_firstdiff_lag6 = lagmatrix(ECM2.Fut18_firstdiff_lag1,5);
ECM2.Fut18_firstdiff_lag7 = lagmatrix(ECM2.Fut18_firstdiff_lag1,6);
ECM2.Fut18_firstdiff_lag8 = lagmatrix(ECM2.Fut18_firstdiff_lag1,7);
%Lag variables spot:
ECM2.Spot_firstdiff_lag1 = lagmatrix(ECM2.Spot_firstdiff,1);
ECM2.Spot_firstdiff_lag2 = lagmatrix(ECM2.Spot_firstdiff,2);
ECM2.Spot_firstdiff_lag3 = lagmatrix(ECM2.Spot_firstdiff,3);
ECM2.Spot_firstdiff_lag4 = lagmatrix(ECM2.Spot_firstdiff,4);
ECM2.Spot_firstdiff_lag5 = lagmatrix(ECM2.Spot_firstdiff,5);
ECM2.Spot_firstdiff_lag6 = lagmatrix(ECM2.Spot_firstdiff,6);
ECM2.Spot_firstdiff_lag7 = lagmatrix(ECM2.Spot_firstdiff,7);
ECM2.Spot_firstdiff_lag8 = lagmatrix(ECM2.Spot_firstdiff,8);


ECM_2 = fitlm(ECM2, ['Spot_firstdiff ~ LaggedResiduals + Fut18_firstdiff_lag1 + Fut18_firstdiff_lag2' ... 
    '+ Fut18_firstdiff_lag3 + Fut18_firstdiff_lag4 + Fut18_firstdiff_lag5 + Fut18_firstdiff_lag6' ...
    '+ Fut18_firstdiff_lag7 + Fut18_firstdiff_lag8 + Spot_firstdiff_lag1 + Spot_firstdiff_lag2' ...
    '+ Spot_firstdiff_lag3 + Spot_firstdiff_lag4 + Spot_firstdiff_lag5 + Spot_firstdiff_lag6' ...
    '+ Spot_firstdiff_lag7 + Spot_firstdiff_lag8'])

% Engle's ARCH test - Spot and Fut18:
ECM2.SqResiduals = (ECM_2.Residuals.Raw).^2;
ECM2.res1 = lagmatrix(ECM2.SqResiduals,1);
ECM2.res2 = lagmatrix(ECM2.SqResiduals,2);
ECM2.res3 = lagmatrix(ECM2.SqResiduals,3);
ECM2.res4 = lagmatrix(ECM2.SqResiduals,4);
ECM2.res5 = lagmatrix(ECM2.SqResiduals,5);
ECM2.res6 = lagmatrix(ECM2.SqResiduals,6);
ECM2.res7 = lagmatrix(ECM2.SqResiduals,7);
ECM2.res8 = lagmatrix(ECM2.SqResiduals,8);
ECM2.res9 = lagmatrix(ECM2.SqResiduals,9);
ECM2.res10 = lagmatrix(ECM2.SqResiduals,10);
ECM2.res11 = lagmatrix(ECM2.SqResiduals,11);
ECM2.res12 = lagmatrix(ECM2.SqResiduals,12);

ECM2 = fillmissing(ECM2, 'constant', 0, 'DataVariables', ...
                    {'res1','res2','res3','res4','res5','res6','res7', ...
                    'res8','res9','res10','res11','res12'});

% Running the auxiliary regression
AuxReg2 = fitlm(ECM2, ['SqResiduals ~ res1 + res2 + res3 + res3 + res4 + res5' ...
    '+ res6 + res7 + res8 + res9 + res10 + res11 + res12'])

% Testing for autocorrelation:
T =  1446;
alpha = [0.1, 0.05, 0.01]
TestStatARCH = (T)*AuxReg2.Rsquared.Ordinary
pVal = 1 - cdf('Chisquare', TestStatARCH, 12)
criticalARCH = chi2inv(1-alpha, 12);
rejection = TestStatARCH > criticalARCH

% ECM-GARCH: Spot and Fut18:
GARCH2 = garch('GARCHlags', 1, 'ARCHlags', 1)
[estMdl, estParamCov, logL] = estimate(GARCH2, ECM2.Spot_firstdiff)
condVar = infer(estMdl, ECM2.Spot_firstdiff);
ECM2.condVol = sqrt(condVar);

ECM_GARCH_2 = fitlm(ECM2, ['Spot_firstdiff ~ LaggedResiduals + Fut18_firstdiff_lag1 + Fut18_firstdiff_lag2' ... 
    '+ Fut18_firstdiff_lag3 + Fut18_firstdiff_lag4 + Fut18_firstdiff_lag5 + Fut18_firstdiff_lag6' ...
    '+ Fut18_firstdiff_lag7 + Fut18_firstdiff_lag8 + Spot_firstdiff_lag1 + Spot_firstdiff_lag2' ...
    '+ Spot_firstdiff_lag3 + Spot_firstdiff_lag4 + Spot_firstdiff_lag5 + Spot_firstdiff_lag6' ...
    '+ Spot_firstdiff_lag7 + Spot_firstdiff_lag8 + condVol'])

figure;
plot(ECM2.Date, ECM2.Spot_firstdiff);
hold on
plot(ECM2.Date, ECM2.condVol);
hold off
legend('Log-differenced spot price', 'Inferred volatility')
box('off')

% Finding residuals between Spot and lagged Fut17:
reg3 = table();
reg3.Date = data.Date;
reg3.Spot = data.Spot;
reg3.LagFut17 = lagmatrix(data.Fut17,1);
reg3(1,:) = [];
reg3(1208:end,:) = [];
model3 = fitlm(reg3, 'Spot ~ LagFut17')
reg3.residuals = model3.Residuals.Raw;
plot(reg3.Date, reg3.residuals)

% Defining ECM variables - Spot and Fut17:
ECM3 = table();
ECM3.Date = reg3.Date;
ECM3.Spot_firstdiff = data.Spot_firstdiff(2:1208);
ECM3.LaggedResiduals = lagmatrix(reg3.residuals,1);
%Lag variables - Fut:
ECM3.Fut17_firstdiff_lag1 = reg3.LagFut17 - lagmatrix(reg3.LagFut17,1);
ECM3.Fut17_firstdiff_lag2 = lagmatrix(ECM3.Fut17_firstdiff_lag1,1);
ECM3.Fut17_firstdiff_lag3 = lagmatrix(ECM3.Fut17_firstdiff_lag1,2);
ECM3.Fut17_firstdiff_lag4 = lagmatrix(ECM3.Fut17_firstdiff_lag1,3);
ECM3.Fut17_firstdiff_lag5 = lagmatrix(ECM3.Fut17_firstdiff_lag1,4);
ECM3.Fut17_firstdiff_lag6 = lagmatrix(ECM3.Fut17_firstdiff_lag1,5);
ECM3.Fut17_firstdiff_lag7 = lagmatrix(ECM3.Fut17_firstdiff_lag1,6);
ECM3.Fut17_firstdiff_lag8 = lagmatrix(ECM3.Fut17_firstdiff_lag1,7);
ECM3.Fut17_firstdiff_lag9 = lagmatrix(ECM3.Fut17_firstdiff_lag1,8);
ECM3.Fut17_firstdiff_lag10 = lagmatrix(ECM3.Fut17_firstdiff_lag1,9);
ECM3.Fut17_firstdiff_lag11 = lagmatrix(ECM3.Fut17_firstdiff_lag1,10);
ECM3.Fut17_firstdiff_lag12 = lagmatrix(ECM3.Fut17_firstdiff_lag1,11);

%Lag variables spot:
ECM3.Spot_firstdiff_lag1 = lagmatrix(ECM3.Spot_firstdiff,1);
ECM3.Spot_firstdiff_lag2 = lagmatrix(ECM3.Spot_firstdiff,2);
ECM3.Spot_firstdiff_lag3 = lagmatrix(ECM3.Spot_firstdiff,3);
ECM3.Spot_firstdiff_lag4 = lagmatrix(ECM3.Spot_firstdiff,4);
ECM3.Spot_firstdiff_lag5 = lagmatrix(ECM3.Spot_firstdiff,5);
ECM3.Spot_firstdiff_lag6 = lagmatrix(ECM3.Spot_firstdiff,6);
ECM3.Spot_firstdiff_lag7 = lagmatrix(ECM3.Spot_firstdiff,7);
ECM3.Spot_firstdiff_lag8 = lagmatrix(ECM3.Spot_firstdiff,8);
ECM3.Spot_firstdiff_lag9 = lagmatrix(ECM3.Spot_firstdiff,9);
ECM3.Spot_firstdiff_lag10 = lagmatrix(ECM3.Spot_firstdiff,10);
ECM3.Spot_firstdiff_lag11 = lagmatrix(ECM3.Spot_firstdiff,11);
ECM3.Spot_firstdiff_lag12 = lagmatrix(ECM3.Spot_firstdiff,12);


ECM_3 = fitlm(ECM3, ['Spot_firstdiff ~ LaggedResiduals + Fut17_firstdiff_lag1 + Fut17_firstdiff_lag2' ... 
    '+ Fut17_firstdiff_lag3 + Fut17_firstdiff_lag4 + Fut17_firstdiff_lag5 + Fut17_firstdiff_lag6' ...
    '+ Fut17_firstdiff_lag7 + Fut17_firstdiff_lag8 + Fut17_firstdiff_lag9 + Fut17_firstdiff_lag10' ...
    '+ Fut17_firstdiff_lag11 + Fut17_firstdiff_lag12 + Spot_firstdiff_lag1 + Spot_firstdiff_lag2' ...
    '+ Spot_firstdiff_lag3 + Spot_firstdiff_lag4 + Spot_firstdiff_lag5 + Spot_firstdiff_lag6' ...
    '+ Spot_firstdiff_lag7 + Spot_firstdiff_lag8 + Spot_firstdiff_lag9 + Spot_firstdiff_lag10' ...
    '+ Spot_firstdiff_lag11 + Spot_firstdiff_lag12'])

% Engle's ARCH Test - Spot and Fut17:
ECM3.SqResiduals = (ECM_3.Residuals.Raw).^2;
ECM3.res1 = lagmatrix(ECM3.SqResiduals,1);
ECM3.res2 = lagmatrix(ECM3.SqResiduals,2);
ECM3.res3 = lagmatrix(ECM3.SqResiduals,3);
ECM3.res4 = lagmatrix(ECM3.SqResiduals,4);
ECM3.res5 = lagmatrix(ECM3.SqResiduals,5);
ECM3.res6 = lagmatrix(ECM3.SqResiduals,6);
ECM3.res7 = lagmatrix(ECM3.SqResiduals,7);
ECM3.res8 = lagmatrix(ECM3.SqResiduals,8);
ECM3.res9 = lagmatrix(ECM3.SqResiduals,9);
ECM3.res10 = lagmatrix(ECM3.SqResiduals,10);
ECM3.res11 = lagmatrix(ECM3.SqResiduals,11);
ECM3.res12 = lagmatrix(ECM3.SqResiduals,12);

ECM3 = fillmissing(ECM3, 'constant', 0, 'DataVariables', ...
                    {'res1','res2','res3','res4','res5','res6','res7', ...
                    'res8','res9','res10','res11','res12'});

% Running the auxiliary regression
AuxReg3 = fitlm(ECM3, ['SqResiduals ~ res1 + res2 + res3 + res3 + res4 + res5' ...
    '+ res6 + res7 + res8 + res9 + res10 + res11 + res12'])

% Testing for autocorrelation:
T =  1195;
alpha = [0.1, 0.05, 0.01]
TestStatARCH = (T)*AuxReg3.Rsquared.Ordinary
pVal = 1 - cdf('Chisquare', TestStatARCH, 12)
criticalARCH = chi2inv(1-alpha, 12);
rejection = TestStatARCH > criticalARCH

% ECM-GARCH: Spot and Fut17:
GARCH3 = garch('GARCHlags', 1, 'ARCHlags', 1)
[estMdl, estParamCov, logL] = estimate(GARCH3, ECM3.Spot_firstdiff)
condVar = infer(estMdl, ECM3.Spot_firstdiff);
ECM3.condVol = sqrt(condVar);

ECM_GARCH_3 = fitlm(ECM3, ['Spot_firstdiff ~ LaggedResiduals + Fut17_firstdiff_lag1 + Fut17_firstdiff_lag2' ... 
    '+ Fut17_firstdiff_lag3 + Fut17_firstdiff_lag4 + Fut17_firstdiff_lag5 + Fut17_firstdiff_lag6' ...
    '+ Fut17_firstdiff_lag7 + Fut17_firstdiff_lag8 + Fut17_firstdiff_lag9 + Fut17_firstdiff_lag10' ...
    '+ Fut17_firstdiff_lag11 + Fut17_firstdiff_lag12 + Spot_firstdiff_lag1 + Spot_firstdiff_lag2' ...
    '+ Spot_firstdiff_lag3 + Spot_firstdiff_lag4 + Spot_firstdiff_lag5 + Spot_firstdiff_lag6' ...
    '+ Spot_firstdiff_lag7 + Spot_firstdiff_lag8 + Spot_firstdiff_lag9 + Spot_firstdiff_lag10' ...
    '+ Spot_firstdiff_lag11 + Spot_firstdiff_lag12 + condVol'])

figure;
plot(ECM3.Date, ECM3.Spot_firstdiff);
hold on
plot(ECM3.Date, ECM3.condVol);
hold off
legend('Log-differenced spot price', 'Inferred volatility')
box('off')
xtickformat('yyyy')

% Finding residuals between Spot and lagged Fut16:
reg4 = table();
reg4.Date = data.Date;
reg4.Spot = data.Spot;
reg4.LagFut16 = lagmatrix(data.Fut16,1);
reg4(1,:) = [];
reg4(984:end,:) = [];
model4 = fitlm(reg4, 'Spot ~ LagFut16')
reg4.residuals = model4.Residuals.Raw;
plot(reg4.Date, reg4.residuals)

% Defining ECM variables - Spot and Fut16:
ECM4 = table();
ECM4.Date = reg4.Date;
ECM4.Spot_firstdiff = data.Spot_firstdiff(2:984);
ECM4.LaggedResiduals = lagmatrix(reg4.residuals,1);
%Lag variables - Fut:
ECM4.Fut16_firstdiff_lag1 = reg4.LagFut16 - lagmatrix(reg4.LagFut16,1);
ECM4.Fut16_firstdiff_lag2 = lagmatrix(ECM4.Fut16_firstdiff_lag1,1);
ECM4.Fut16_firstdiff_lag3 = lagmatrix(ECM4.Fut16_firstdiff_lag1,2);
ECM4.Fut16_firstdiff_lag4 = lagmatrix(ECM4.Fut16_firstdiff_lag1,3);
ECM4.Fut16_firstdiff_lag5 = lagmatrix(ECM4.Fut16_firstdiff_lag1,4);
ECM4.Fut16_firstdiff_lag6 = lagmatrix(ECM4.Fut16_firstdiff_lag1,5);
ECM4.Fut16_firstdiff_lag7 = lagmatrix(ECM4.Fut16_firstdiff_lag1,6);
%Lag variables spot:
ECM4.Spot_firstdiff_lag1 = lagmatrix(ECM4.Spot_firstdiff,1);
ECM4.Spot_firstdiff_lag2 = lagmatrix(ECM4.Spot_firstdiff,2);
ECM4.Spot_firstdiff_lag3 = lagmatrix(ECM4.Spot_firstdiff,3);
ECM4.Spot_firstdiff_lag4 = lagmatrix(ECM4.Spot_firstdiff,4);
ECM4.Spot_firstdiff_lag5 = lagmatrix(ECM4.Spot_firstdiff,5);
ECM4.Spot_firstdiff_lag6 = lagmatrix(ECM4.Spot_firstdiff,6);
ECM4.Spot_firstdiff_lag7 = lagmatrix(ECM4.Spot_firstdiff,7);

ECM_4 = fitlm(ECM4, ['Spot_firstdiff ~ LaggedResiduals + Fut16_firstdiff_lag1 + Fut16_firstdiff_lag2' ... 
    '+ Fut16_firstdiff_lag3 + Fut16_firstdiff_lag4 + Fut16_firstdiff_lag5 + Fut16_firstdiff_lag6' ...
    '+ Fut16_firstdiff_lag7 + Spot_firstdiff_lag1 + Spot_firstdiff_lag2' ...
    '+ Spot_firstdiff_lag3 + Spot_firstdiff_lag4 + Spot_firstdiff_lag5 + Spot_firstdiff_lag6' ...
    '+ Spot_firstdiff_lag7'])

% Engle's ARCH Test - Spot and Fut16:
ECM4.SqResiduals = (ECM_4.Residuals.Raw).^2;
ECM4.res1 = lagmatrix(ECM4.SqResiduals,1);
ECM4.res2 = lagmatrix(ECM4.SqResiduals,2);
ECM4.res3 = lagmatrix(ECM4.SqResiduals,3);
ECM4.res4 = lagmatrix(ECM4.SqResiduals,4);
ECM4.res5 = lagmatrix(ECM4.SqResiduals,5);
ECM4.res6 = lagmatrix(ECM4.SqResiduals,6);
ECM4.res7 = lagmatrix(ECM4.SqResiduals,7);
ECM4.res8 = lagmatrix(ECM4.SqResiduals,8);
ECM4.res9 = lagmatrix(ECM4.SqResiduals,9);
ECM4.res10 = lagmatrix(ECM4.SqResiduals,10);
ECM4.res11 = lagmatrix(ECM4.SqResiduals,11);
ECM4.res12 = lagmatrix(ECM4.SqResiduals,12);

ECM4 = fillmissing(ECM4, 'constant', 0, 'DataVariables', ...
                    {'res1','res2','res3','res4','res5','res6','res7', ...
                    'res8','res9','res10','res11','res12'});

% Running the auxiliary regression
AuxReg4 = fitlm(ECM4, ['SqResiduals ~ res1 + res2 + res3 + res3 + res4 + res5' ...
    '+ res6 + res7 + res8 + res9 + res10 + res11 + res12'])

% Testing for autocorrelation:
T =  976;
alpha = [0.1, 0.05, 0.01]
TestStatARCH = (T)*AuxReg4.Rsquared.Ordinary
pVal = 1 - cdf('Chisquare', TestStatARCH, 12)
criticalARCH = chi2inv(1-alpha, 12);
rejection = TestStatARCH > criticalARCH

% ECM-GARCH: Spot and Fut16:
GARCH4 = garch('GARCHlags', 1, 'ARCHlags', 1)
[estMdl, estParamCov, logL] = estimate(GARCH4, ECM4.Spot_firstdiff)
condVar = infer(estMdl, ECM4.Spot_firstdiff);
ECM4.condVol = sqrt(condVar);

ECM_GARCH_4 = fitlm(ECM4, ['Spot_firstdiff ~ LaggedResiduals + Fut16_firstdiff_lag1 + Fut16_firstdiff_lag2' ... 
    '+ Fut16_firstdiff_lag3 + Fut16_firstdiff_lag4 + Fut16_firstdiff_lag5 + Fut16_firstdiff_lag6' ...
    '+ Fut16_firstdiff_lag7 + Spot_firstdiff_lag1 + Spot_firstdiff_lag2' ...
    '+ Spot_firstdiff_lag3 + Spot_firstdiff_lag4 + Spot_firstdiff_lag5 + Spot_firstdiff_lag6' ...
    '+ Spot_firstdiff_lag7 + condVol'])

figure;
plot(ECM4.Date, ECM4.Spot_firstdiff);
hold on
plot(ECM4.Date, ECM4.condVol)
hold off
legend('Log-differenced spot price', 'Inferred volatility')
box('off')
xtickformat('yyyy')

% Finding residuals between Spot and lagged Fut15:
reg5 = table();
reg5.Date = data.Date;
reg5.Spot = data.Spot;
reg5.LagFut15 = lagmatrix(data.Fut15,1);
reg5(1,:) = [];
reg5(730:end,:) = [];
model5 = fitlm(reg5, 'Spot ~ LagFut15')
reg5.residuals = model5.Residuals.Raw;
plot(reg5.Date, reg5.residuals)

% Defining ECM variables - Spot and Fut15:
ECM5 = table();
ECM5.Date = reg5.Date;
ECM5.Spot_firstdiff = data.Spot_firstdiff(2:730);
ECM5.LaggedResiduals = lagmatrix(reg5.residuals,1);
%Lag variables - Fut:
ECM5.Fut15_firstdiff_lag1 = reg5.LagFut15 - lagmatrix(reg5.LagFut15,1);
ECM5.Fut15_firstdiff_lag2 = lagmatrix(ECM5.Fut15_firstdiff_lag1,1);
ECM5.Fut15_firstdiff_lag3 = lagmatrix(ECM5.Fut15_firstdiff_lag1,2);
ECM5.Fut15_firstdiff_lag4 = lagmatrix(ECM5.Fut15_firstdiff_lag1,3);
ECM5.Fut15_firstdiff_lag5 = lagmatrix(ECM5.Fut15_firstdiff_lag1,4);
ECM5.Fut15_firstdiff_lag6 = lagmatrix(ECM5.Fut15_firstdiff_lag1,5);
ECM5.Fut15_firstdiff_lag7 = lagmatrix(ECM5.Fut15_firstdiff_lag1,6);
ECM5.Fut15_firstdiff_lag8 = lagmatrix(ECM5.Fut15_firstdiff_lag1,7);
ECM5.Fut15_firstdiff_lag9 = lagmatrix(ECM5.Fut15_firstdiff_lag1,8);
ECM5.Fut15_firstdiff_lag10 = lagmatrix(ECM5.Fut15_firstdiff_lag1,9);

%Lag variables spot:
ECM5.Spot_firstdiff_lag1 = lagmatrix(ECM5.Spot_firstdiff,1);
ECM5.Spot_firstdiff_lag2 = lagmatrix(ECM5.Spot_firstdiff,2);
ECM5.Spot_firstdiff_lag3 = lagmatrix(ECM5.Spot_firstdiff,3);
ECM5.Spot_firstdiff_lag4 = lagmatrix(ECM5.Spot_firstdiff,4);
ECM5.Spot_firstdiff_lag5 = lagmatrix(ECM5.Spot_firstdiff,5);
ECM5.Spot_firstdiff_lag6 = lagmatrix(ECM5.Spot_firstdiff,6);
ECM5.Spot_firstdiff_lag7 = lagmatrix(ECM5.Spot_firstdiff,7);
ECM5.Spot_firstdiff_lag8 = lagmatrix(ECM5.Spot_firstdiff,8);
ECM5.Spot_firstdiff_lag9 = lagmatrix(ECM5.Spot_firstdiff,9);
ECM5.Spot_firstdiff_lag10 = lagmatrix(ECM5.Spot_firstdiff,10);



ECM_5 = fitlm(ECM5, ['Spot_firstdiff ~ LaggedResiduals + Fut15_firstdiff_lag1 + Fut15_firstdiff_lag2' ... 
    '+ Fut15_firstdiff_lag3 + Fut15_firstdiff_lag4 + Fut15_firstdiff_lag5 + Fut15_firstdiff_lag6' ...
    '+ Fut15_firstdiff_lag7 + Fut15_firstdiff_lag8 + Fut15_firstdiff_lag9 + Fut15_firstdiff_lag10' ...
    '+ Spot_firstdiff_lag1 + Spot_firstdiff_lag2' ...
    '+ Spot_firstdiff_lag3 + Spot_firstdiff_lag4 + Spot_firstdiff_lag5 + Spot_firstdiff_lag6' ...
    '+ Spot_firstdiff_lag7 + Spot_firstdiff_lag8 + Spot_firstdiff_lag9 + Spot_firstdiff_lag10'])

% Engle's ARCH Test - Spot and Fut15:
ECM5.SqResiduals = (ECM_5.Residuals.Raw).^2;
ECM5.res1 = lagmatrix(ECM5.SqResiduals,1);
ECM5.res2 = lagmatrix(ECM5.SqResiduals,2);
ECM5.res3 = lagmatrix(ECM5.SqResiduals,3);
ECM5.res4 = lagmatrix(ECM5.SqResiduals,4);
ECM5.res5 = lagmatrix(ECM5.SqResiduals,5);
ECM5.res6 = lagmatrix(ECM5.SqResiduals,6);
ECM5.res7 = lagmatrix(ECM5.SqResiduals,7);
ECM5.res8 = lagmatrix(ECM5.SqResiduals,8);
ECM5.res9 = lagmatrix(ECM5.SqResiduals,9);
ECM5.res10 = lagmatrix(ECM5.SqResiduals,10);
ECM5.res11 = lagmatrix(ECM5.SqResiduals,11);
ECM5.res12 = lagmatrix(ECM5.SqResiduals,12);

ECM5 = fillmissing(ECM5, 'constant', 0, 'DataVariables', ...
                    {'res1','res2','res3','res4','res5','res6','res7', ...
                    'res8','res9','res10','res11','res12'});

% Running the auxiliary regression
AuxReg5 = fitlm(ECM5, ['SqResiduals ~ res1 + res2 + res3 + res3 + res4 + res5' ...
    '+ res6 + res7 + res8 + res9 + res10 + res11 + res12'])

% Testing for ARCH effects:
T =  719;
alpha = [0.1, 0.05, 0.01]
TestStatARCH = (T)*AuxReg5.Rsquared.Ordinary
pVal = 1 - cdf('Chisquare', TestStatARCH, 12)
criticalARCH = chi2inv(1-alpha, 12);
rejection = TestStatARCH > criticalARCH

% ECM-GARCH: Spot and Fut15:
GARCH5 = garch('GARCHlags', 1, 'ARCHlags', 1)
[estMdl, estParamCov, logL] = estimate(GARCH5, ECM5.Spot_firstdiff)
condVar = infer(estMdl, ECM5.Spot_firstdiff);
ECM5.condVol = sqrt(condVar);

ECM_GARCH_5 = fitlm(ECM5, ['Spot_firstdiff ~ LaggedResiduals + Fut15_firstdiff_lag1 + Fut15_firstdiff_lag2' ... 
    '+ Fut15_firstdiff_lag3 + Fut15_firstdiff_lag4 + Fut15_firstdiff_lag5 + Fut15_firstdiff_lag6' ...
    '+ Fut15_firstdiff_lag7 + Fut15_firstdiff_lag8 + Fut15_firstdiff_lag9 + Fut15_firstdiff_lag10' ...
    '+ Spot_firstdiff_lag1 + Spot_firstdiff_lag2 + Spot_firstdiff_lag3 + Spot_firstdiff_lag4' ...
    '+ Spot_firstdiff_lag5 + Spot_firstdiff_lag6+ Spot_firstdiff_lag7 + Spot_firstdiff_lag8' ...
    '+ Spot_firstdiff_lag9 + Spot_firstdiff_lag10 + condVol'])

figure;
plot(ECM5.Date, ECM5.Spot_firstdiff);
hold on
plot(ECM5.Date, ECM5.condVol);
hold off
legend('Log-differenced spot price', 'Inferred volatility')
box('off')
xtickformat('yyyy')

% Finding residuals between Spot and lagged Fut14:
reg6 = table();
reg6.Date = data.Date;
reg6.Spot = data.Spot;
reg6.LagFut14 = lagmatrix(data.Fut14,1);
reg6(1,:) = [];
reg6(486:end,:) = [];
model6 = fitlm(reg6, 'Spot ~ LagFut14')
reg6.residuals = model6.Residuals.Raw;
plot(reg6.Date, reg6.residuals)

% Defining ECM variables - Spot and Fut14:
ECM6 = table();
ECM6.Date = reg6.Date;
ECM6.Spot_firstdiff = data.Spot_firstdiff(2:486);
ECM6.LaggedResiduals = lagmatrix(reg6.residuals,1);
%Lag variables - Fut:
ECM6.Fut14_firstdiff_lag1 = reg6.LagFut14 - lagmatrix(reg6.LagFut14,1);
ECM6.Fut14_firstdiff_lag2 = lagmatrix(ECM6.Fut14_firstdiff_lag1,1);
ECM6.Fut14_firstdiff_lag3 = lagmatrix(ECM6.Fut14_firstdiff_lag1,2);
ECM6.Fut14_firstdiff_lag4 = lagmatrix(ECM6.Fut14_firstdiff_lag1,3);
ECM6.Fut14_firstdiff_lag5 = lagmatrix(ECM6.Fut14_firstdiff_lag1,4);
ECM6.Fut14_firstdiff_lag6 = lagmatrix(ECM6.Fut14_firstdiff_lag1,5);
ECM6.Fut14_firstdiff_lag7 = lagmatrix(ECM6.Fut14_firstdiff_lag1,6);

%Lag variables spot:
ECM6.Spot_firstdiff_lag1 = lagmatrix(ECM6.Spot_firstdiff,1);
ECM6.Spot_firstdiff_lag2 = lagmatrix(ECM6.Spot_firstdiff,2);
ECM6.Spot_firstdiff_lag3 = lagmatrix(ECM6.Spot_firstdiff,3);
ECM6.Spot_firstdiff_lag4 = lagmatrix(ECM6.Spot_firstdiff,4);
ECM6.Spot_firstdiff_lag5 = lagmatrix(ECM6.Spot_firstdiff,5);
ECM6.Spot_firstdiff_lag6 = lagmatrix(ECM6.Spot_firstdiff,6);
ECM6.Spot_firstdiff_lag7 = lagmatrix(ECM6.Spot_firstdiff,7);

ECM_6 = fitlm(ECM6, ['Spot_firstdiff ~ LaggedResiduals + Fut14_firstdiff_lag1 + Fut14_firstdiff_lag2' ... 
    '+ Fut14_firstdiff_lag3 + Fut14_firstdiff_lag4 + Fut14_firstdiff_lag5 + Fut14_firstdiff_lag6' ...
    '+ Fut14_firstdiff_lag7 + Spot_firstdiff_lag1 + Spot_firstdiff_lag2' ...
    '+ Spot_firstdiff_lag3 + Spot_firstdiff_lag4 + Spot_firstdiff_lag5 + Spot_firstdiff_lag6' ...
    '+ Spot_firstdiff_lag7'])

% Engle's ARCH Test - Spot and Fut14:
ECM6.SqResiduals = (ECM_6.Residuals.Raw).^2;
ECM6.res1 = lagmatrix(ECM6.SqResiduals,1);
ECM6.res2 = lagmatrix(ECM6.SqResiduals,2);
ECM6.res3 = lagmatrix(ECM6.SqResiduals,3);
ECM6.res4 = lagmatrix(ECM6.SqResiduals,4);
ECM6.res5 = lagmatrix(ECM6.SqResiduals,5);
ECM6.res6 = lagmatrix(ECM6.SqResiduals,6);
ECM6.res7 = lagmatrix(ECM6.SqResiduals,7);
ECM6.res8 = lagmatrix(ECM6.SqResiduals,8);
ECM6.res9 = lagmatrix(ECM6.SqResiduals,9);
ECM6.res10 = lagmatrix(ECM6.SqResiduals,10);
ECM6.res11 = lagmatrix(ECM6.SqResiduals,11);
ECM6.res12 = lagmatrix(ECM6.SqResiduals,12);

ECM6 = fillmissing(ECM6, 'constant', 0, 'DataVariables', ...
                    {'res1','res2','res3','res4','res5','res6','res7', ...
                    'res8','res9','res10','res11','res12'});

% Running the auxiliary regression
AuxReg6 = fitlm(ECM6, ['SqResiduals ~ res1 + res2 + res3 + res3 + res4 + res5' ...
    '+ res6 + res7 + res8 + res9 + res10 + res11 + res12'])

% Testing for ARCH effects:
T =  478;
alpha = [0.1, 0.05, 0.01]
TestStatARCH = (T)*AuxReg6.Rsquared.Ordinary
pVal = 1 - cdf('Chisquare', TestStatARCH, 12)
criticalARCH = chi2inv(1-alpha, 12);
rejection = TestStatARCH > criticalARCH

% ECM-GARCH: Spot and Fut14:
GARCH6 = garch('GARCHlags', 1, 'ARCHlags', 1)
[estMdl, estParamCov, logL] = estimate(GARCH6, ECM6.Spot_firstdiff)
condVar = infer(estMdl, ECM6.Spot_firstdiff);
ECM6.condVol = sqrt(condVar);

ECM_GARCH_6 = fitlm(ECM6, ['Spot_firstdiff ~ LaggedResiduals + Fut14_firstdiff_lag1 + Fut14_firstdiff_lag2' ... 
    '+ Fut14_firstdiff_lag3 + Fut14_firstdiff_lag4 + Fut14_firstdiff_lag5 + Fut14_firstdiff_lag6' ...
    '+ Fut14_firstdiff_lag7 + Spot_firstdiff_lag1 + Spot_firstdiff_lag2' ...
    '+ Spot_firstdiff_lag3 + Spot_firstdiff_lag4 + Spot_firstdiff_lag5 + Spot_firstdiff_lag6' ...
    '+ Spot_firstdiff_lag7 + condVol'])

figure;
plot(ECM6.Date, ECM6.Spot_firstdiff);
hold on
plot(ECM6.Date, ECM6.condVol)
hold off
legend('Log-differenced spot price', 'Inferred volatility')
box('off');
xtickformat('yyyy')

%% Part 5 - Predictions by ECM and ECM-GARCH (re-estimation of all models):
clear;
clc;

% Importing data and creating variables for sample matching after estimation:
data = readtable('Spot+DecFutures.xlsx');
data=data(~any(ismissing(data.Spot),2),:);
data(584:585,:) = [];
data = table2timetable(data);
data.Spot = log(data.Spot);
data.Fut19 = log(data.Fut19);
data.Fut18 = log(data.Fut18);
data.Fut17 = log(data.Fut17);
data.Fut16 = log(data.Fut16);
data.Fut15 = log(data.Fut15);
data.Fut14 = log(data.Fut14);

% Creating first differences and lags:
data.Spot_firstdiff = data.Spot - lagmatrix(data.Spot,1);
data.Spot_firstdiff_lag1 = lagmatrix(data.Spot_firstdiff,1);
data.Spot_firstdiff_lag2 = lagmatrix(data.Spot_firstdiff,2);
data.Spot_firstdiff_lag3 = lagmatrix(data.Spot_firstdiff,3);
data.Spot_firstdiff_lag4 = lagmatrix(data.Spot_firstdiff,4);
data.Spot_firstdiff_lag5 = lagmatrix(data.Spot_firstdiff,5);
data.Spot_firstdiff_lag6 = lagmatrix(data.Spot_firstdiff,6);
data.Spot_firstdiff_lag7 = lagmatrix(data.Spot_firstdiff,7);
data.Spot_firstdiff_lag8 = lagmatrix(data.Spot_firstdiff,8);
data.Spot_firstdiff_lag9 = lagmatrix(data.Spot_firstdiff,9);
data.Spot_firstdiff_lag10 = lagmatrix(data.Spot_firstdiff,10);
data.Spot_firstdiff_lag11 = lagmatrix(data.Spot_firstdiff,11);
data.Spot_firstdiff_lag12 = lagmatrix(data.Spot_firstdiff,12);

data.Fut19_firstdiff = data.Fut19 - lagmatrix(data.Fut19,1);
data.Fut19_firstdiff_lag1 = lagmatrix(data.Fut19_firstdiff,1);
data.Fut19_firstdiff_lag2 = lagmatrix(data.Fut19_firstdiff,2);
data.Fut19_firstdiff_lag3 = lagmatrix(data.Fut19_firstdiff,3);
data.Fut19_firstdiff_lag4 = lagmatrix(data.Fut19_firstdiff,4);
data.Fut19_firstdiff_lag5 = lagmatrix(data.Fut19_firstdiff,5);
data.Fut19_firstdiff_lag6 = lagmatrix(data.Fut19_firstdiff,6);
data.Fut19_firstdiff_lag7 = lagmatrix(data.Fut19_firstdiff,7);
data.Fut19_firstdiff_lag8 = lagmatrix(data.Fut19_firstdiff,8);
data.Fut19_firstdiff_lag9 = lagmatrix(data.Fut19_firstdiff,9);
data.Fut19_firstdiff_lag10 = lagmatrix(data.Fut19_firstdiff,10);
data.Fut19_firstdiff_lag11 = lagmatrix(data.Fut19_firstdiff,11);
data.Fut19_firstdiff_lag12 = lagmatrix(data.Fut19_firstdiff,12);

data.Fut18_firstdiff = data.Fut18 - lagmatrix(data.Fut18,1);
data.Fut18_firstdiff_lag1 = lagmatrix(data.Fut18_firstdiff,1);
data.Fut18_firstdiff_lag2 = lagmatrix(data.Fut18_firstdiff,2);
data.Fut18_firstdiff_lag3 = lagmatrix(data.Fut18_firstdiff,3);
data.Fut18_firstdiff_lag4 = lagmatrix(data.Fut18_firstdiff,4);
data.Fut18_firstdiff_lag5 = lagmatrix(data.Fut18_firstdiff,5);
data.Fut18_firstdiff_lag6 = lagmatrix(data.Fut18_firstdiff,6);
data.Fut18_firstdiff_lag7 = lagmatrix(data.Fut18_firstdiff,7);
data.Fut18_firstdiff_lag8 = lagmatrix(data.Fut18_firstdiff,8);

data.Fut17_firstdiff = data.Fut17 - lagmatrix(data.Fut17,1);
data.Fut17_firstdiff_lag1 = lagmatrix(data.Fut17_firstdiff,1);
data.Fut17_firstdiff_lag2 = lagmatrix(data.Fut17_firstdiff,2);
data.Fut17_firstdiff_lag3 = lagmatrix(data.Fut17_firstdiff,3);
data.Fut17_firstdiff_lag4 = lagmatrix(data.Fut17_firstdiff,4);
data.Fut17_firstdiff_lag5 = lagmatrix(data.Fut17_firstdiff,5);
data.Fut17_firstdiff_lag6 = lagmatrix(data.Fut17_firstdiff,6);
data.Fut17_firstdiff_lag7 = lagmatrix(data.Fut17_firstdiff,7);
data.Fut17_firstdiff_lag8 = lagmatrix(data.Fut17_firstdiff,8);
data.Fut17_firstdiff_lag9 = lagmatrix(data.Fut17_firstdiff,9);
data.Fut17_firstdiff_lag10 = lagmatrix(data.Fut17_firstdiff,10);
data.Fut17_firstdiff_lag11 = lagmatrix(data.Fut17_firstdiff,11);
data.Fut17_firstdiff_lag12 = lagmatrix(data.Fut17_firstdiff,12);

data.Fut16_firstdiff = data.Fut16 - lagmatrix(data.Fut16,1);
data.Fut16_firstdiff_lag1 = lagmatrix(data.Fut16_firstdiff,1);
data.Fut16_firstdiff_lag2 = lagmatrix(data.Fut16_firstdiff,2);
data.Fut16_firstdiff_lag3 = lagmatrix(data.Fut16_firstdiff,3);
data.Fut16_firstdiff_lag4 = lagmatrix(data.Fut16_firstdiff,4);
data.Fut16_firstdiff_lag5 = lagmatrix(data.Fut16_firstdiff,5);
data.Fut16_firstdiff_lag6 = lagmatrix(data.Fut16_firstdiff,6);
data.Fut16_firstdiff_lag7 = lagmatrix(data.Fut16_firstdiff,7);

data.Fut15_firstdiff = data.Fut15 - lagmatrix(data.Fut15,1);
data.Fut15_firstdiff_lag1 = lagmatrix(data.Fut15_firstdiff,1);
data.Fut15_firstdiff_lag2 = lagmatrix(data.Fut15_firstdiff,2);
data.Fut15_firstdiff_lag3 = lagmatrix(data.Fut15_firstdiff,3);
data.Fut15_firstdiff_lag4 = lagmatrix(data.Fut15_firstdiff,4);
data.Fut15_firstdiff_lag5 = lagmatrix(data.Fut15_firstdiff,5);
data.Fut15_firstdiff_lag6 = lagmatrix(data.Fut15_firstdiff,6);
data.Fut15_firstdiff_lag7 = lagmatrix(data.Fut15_firstdiff,7);
data.Fut15_firstdiff_lag8 = lagmatrix(data.Fut15_firstdiff,8);
data.Fut15_firstdiff_lag9 = lagmatrix(data.Fut15_firstdiff,9);
data.Fut15_firstdiff_lag10 = lagmatrix(data.Fut15_firstdiff,10);

data.Fut14_firstdiff = data.Fut14 - lagmatrix(data.Fut14,1);
data.Fut14_firstdiff_lag1 = lagmatrix(data.Fut14_firstdiff,1);
data.Fut14_firstdiff_lag2 = lagmatrix(data.Fut14_firstdiff,2);
data.Fut14_firstdiff_lag3 = lagmatrix(data.Fut14_firstdiff,3);
data.Fut14_firstdiff_lag4 = lagmatrix(data.Fut14_firstdiff,4);
data.Fut14_firstdiff_lag5 = lagmatrix(data.Fut14_firstdiff,5);
data.Fut14_firstdiff_lag6 = lagmatrix(data.Fut14_firstdiff,6);
data.Fut14_firstdiff_lag7 = lagmatrix(data.Fut14_firstdiff,7);

% Finding residuals between Spot and lagged Fut19:
reg1 = table();
reg1.Date = data.Date(1:1460);
reg1.Spot = data.Spot(1:1460);
reg1.LagFut19 = lagmatrix(data.Fut19(1:1460),1);
reg1(1,:) = [];
model1 = fitlm(reg1, 'Spot ~ LagFut19')
reg1.residuals1 = model1.Residuals.Raw;

% Defining ECM variables - Spot and Fut19:
ECM1 = table();
ECM1.Date = reg1.Date;
ECM1.Spot_firstdiff = data.Spot_firstdiff(2:1460);
ECM1.LaggedResiduals = lagmatrix(reg1.residuals1,1);
%Lag variables - Fut:
ECM1.Fut19_firstdiff_lag1 = reg1.LagFut19 - lagmatrix(reg1.LagFut19,1);
ECM1.Fut19_firstdiff_lag2 = lagmatrix(ECM1.Fut19_firstdiff_lag1,1);
ECM1.Fut19_firstdiff_lag3 = lagmatrix(ECM1.Fut19_firstdiff_lag1,2);
ECM1.Fut19_firstdiff_lag4 = lagmatrix(ECM1.Fut19_firstdiff_lag1,3);
ECM1.Fut19_firstdiff_lag5 = lagmatrix(ECM1.Fut19_firstdiff_lag1,4);
ECM1.Fut19_firstdiff_lag6 = lagmatrix(ECM1.Fut19_firstdiff_lag1,5);
ECM1.Fut19_firstdiff_lag7 = lagmatrix(ECM1.Fut19_firstdiff_lag1,6);
ECM1.Fut19_firstdiff_lag8 = lagmatrix(ECM1.Fut19_firstdiff_lag1,7);
ECM1.Fut19_firstdiff_lag9 = lagmatrix(ECM1.Fut19_firstdiff_lag1,8);
ECM1.Fut19_firstdiff_lag10 = lagmatrix(ECM1.Fut19_firstdiff_lag1,9);
ECM1.Fut19_firstdiff_lag11 = lagmatrix(ECM1.Fut19_firstdiff_lag1,10);
ECM1.Fut19_firstdiff_lag12 = lagmatrix(ECM1.Fut19_firstdiff_lag1,11);
%Lag variables spot:
ECM1.Spot_firstdiff_lag1 = lagmatrix(ECM1.Spot_firstdiff,1);
ECM1.Spot_firstdiff_lag2 = lagmatrix(ECM1.Spot_firstdiff,2);
ECM1.Spot_firstdiff_lag3 = lagmatrix(ECM1.Spot_firstdiff,3);
ECM1.Spot_firstdiff_lag4 = lagmatrix(ECM1.Spot_firstdiff,4);
ECM1.Spot_firstdiff_lag5 = lagmatrix(ECM1.Spot_firstdiff,5);
ECM1.Spot_firstdiff_lag6 = lagmatrix(ECM1.Spot_firstdiff,6);
ECM1.Spot_firstdiff_lag7 = lagmatrix(ECM1.Spot_firstdiff,7);
ECM1.Spot_firstdiff_lag8 = lagmatrix(ECM1.Spot_firstdiff,8);
ECM1.Spot_firstdiff_lag9 = lagmatrix(ECM1.Spot_firstdiff,9);
ECM1.Spot_firstdiff_lag10 = lagmatrix(ECM1.Spot_firstdiff,10);
ECM1.Spot_firstdiff_lag11 = lagmatrix(ECM1.Spot_firstdiff,11);
ECM1.Spot_firstdiff_lag12 = lagmatrix(ECM1.Spot_firstdiff,12);

ECM_1 = fitlm(ECM1, ['Spot_firstdiff ~ LaggedResiduals + Fut19_firstdiff_lag1 + Fut19_firstdiff_lag2' ... 
    '+ Fut19_firstdiff_lag3 + Fut19_firstdiff_lag4 + Fut19_firstdiff_lag5 + Fut19_firstdiff_lag6' ...
    '+ Fut19_firstdiff_lag7 + Fut19_firstdiff_lag8 + Fut19_firstdiff_lag9 + Fut19_firstdiff_lag10' ...
    '+ Fut19_firstdiff_lag11 + Fut19_firstdiff_lag12 + Spot_firstdiff_lag1 + Spot_firstdiff_lag2' ...
    '+ Spot_firstdiff_lag3 + Spot_firstdiff_lag4 + Spot_firstdiff_lag5 + Spot_firstdiff_lag6' ...
    '+ Spot_firstdiff_lag7 + Spot_firstdiff_lag8 + Spot_firstdiff_lag9 + Spot_firstdiff_lag10' ...
    '+ Spot_firstdiff_lag11 + Spot_firstdiff_lag12'])

residuals = table();
residuals.Date = data.Date;
residuals.LaggedResiduals = zeros(1709,1);

for i = 1460:1709
   reg = table();
   reg.Date = data.Date(1:i);
   reg.Spot = data.Spot(1:i);
   reg.LagFut19 = lagmatrix(data.Fut19(1:i),1);
   reg(1,:) = [];
   model = fitlm(reg, 'Spot ~ LagFut19');
   reg.residuals = model.Residuals.Raw; 
   residuals.LaggedResiduals(i+1) = reg.residuals(end);
end

residuals(1:1460,:) =[]; 
residuals(250,:) = [];
    
Predicted_SpotFD_ECM1 = [ECM_1.Coefficients.Estimate(1) + ECM_1.Coefficients.Estimate(3)*data.Fut19_firstdiff_lag1(1461:end) ...
    + ECM_1.Coefficients.Estimate(4)*data.Fut19_firstdiff_lag2(1461:end) + ECM_1.Coefficients.Estimate(5)*data.Fut19_firstdiff_lag3(1461:end) ...
    + ECM_1.Coefficients.Estimate(6)*data.Fut19_firstdiff_lag4(1461:end) + ECM_1.Coefficients.Estimate(7)*data.Fut19_firstdiff_lag5(1461:end) ...
    + ECM_1.Coefficients.Estimate(8)*data.Fut19_firstdiff_lag6(1461:end) + ECM_1.Coefficients.Estimate(9)*data.Fut19_firstdiff_lag7(1461:end) ...
    + ECM_1.Coefficients.Estimate(10)*data.Fut19_firstdiff_lag8(1461:end) + ECM_1.Coefficients.Estimate(11)*data.Fut19_firstdiff_lag9(1461:end) ...
    + ECM_1.Coefficients.Estimate(12)*data.Fut19_firstdiff_lag10(1461:end) + ECM_1.Coefficients.Estimate(13)*data.Fut19_firstdiff_lag11(1461:end) ...
    + ECM_1.Coefficients.Estimate(14)*data.Fut19_firstdiff_lag12(1461:end) + ECM_1.Coefficients.Estimate(15)*data.Spot_firstdiff_lag1(1461:end) ...
    + ECM_1.Coefficients.Estimate(16)*data.Spot_firstdiff_lag2(1461:end) + ECM_1.Coefficients.Estimate(17)*data.Spot_firstdiff_lag3(1461:end) ...
    + ECM_1.Coefficients.Estimate(18)*data.Spot_firstdiff_lag4(1461:end) + ECM_1.Coefficients.Estimate(19)*data.Spot_firstdiff_lag5(1461:end) ...
    + ECM_1.Coefficients.Estimate(20)*data.Spot_firstdiff_lag6(1461:end) + ECM_1.Coefficients.Estimate(21)*data.Spot_firstdiff_lag7(1461:end) ...
    + ECM_1.Coefficients.Estimate(22)*data.Spot_firstdiff_lag8(1461:end) + ECM_1.Coefficients.Estimate(23)*data.Spot_firstdiff_lag9(1461:end) ...
    + ECM_1.Coefficients.Estimate(24)*data.Spot_firstdiff_lag10(1461:end) + ECM_1.Coefficients.Estimate(25)*data.Spot_firstdiff_lag11(1461:end) ...
    + ECM_1.Coefficients.Estimate(26)*data.Spot_firstdiff_lag12(1461:end) + ECM_1.Coefficients.Estimate(2)*residuals.LaggedResiduals];

Predictions_ECM1 = table();
Predictions_ECM1.Date = data.Date(1460:end);
Predictions_ECM1.Spot = data.Spot(1460:end);
Predictions_ECM1.PredictedSpotFD = zeros(250,1);
Predictions_ECM1.PredictedSpotFD(2:end) = Predicted_SpotFD_ECM1;
Predictions_ECM1.PredictedSpot = lagmatrix(Predictions_ECM1.Spot, 1) + Predictions_ECM1.PredictedSpotFD;
Predictions_ECM1(1,:) = [];
Predictions_ECM1.Spot = exp(Predictions_ECM1.Spot);
Predictions_ECM1.PredictedSpot = exp(Predictions_ECM1.PredictedSpot);
Predictions_ECM1(241:end,:) = [];

MSE_ECM1 = (sum(Predictions_ECM1.PredictedSpot - Predictions_ECM1.Spot)^2)/(240);
RMSE_ECM1 = sqrt(MSE_ECM1)

ChangePred_ECM1 = Predictions_ECM1.PredictedSpot - lagmatrix(Predictions_ECM1.PredictedSpot, 1);
ChangeSpot_ECM1 = Predictions_ECM1.Spot - lagmatrix(Predictions_ECM1.Spot, 1);

SignPred_ECM1 = table();
for i = 1:height(Predictions_ECM1)
    if ChangePred_ECM1(i) >= 0 & ChangeSpot_ECM1(i) >= 0 
        SignPred_ECM1.True(i) = 1;
    else
        SignPred_ECM1.True(i) = 0;
    end
end

for i = 1:height(Predictions_ECM1)
    if ChangePred_ECM1(i) <= 0 & ChangeSpot_ECM1(i) <= 0 
        SignPred_ECM1.True1(i) = 1;
    else
        SignPred_ECM1.True1(i) = 0;
    end
end
PercentCorrect_ECM1 = (sum(SignPred_ECM1.True) + sum(SignPred_ECM1.True1)) / (height(Predictions_ECM1));

% ECM-GARCH: Spot and Fut19:
GARCH1 = garch('GARCHlags', 1, 'ARCHlags', 1)
[estMdl, estParamCov, logL] = estimate(GARCH1, ECM1.Spot_firstdiff)
condVar = infer(estMdl, ECM1.Spot_firstdiff);
ECM1.condVol = sqrt(condVar);

ECM_GARCH_1 = fitlm(ECM1, ['Spot_firstdiff ~ LaggedResiduals + Fut19_firstdiff_lag1 + Fut19_firstdiff_lag2' ... 
    '+ Fut19_firstdiff_lag3 + Fut19_firstdiff_lag4 + Fut19_firstdiff_lag5 + Fut19_firstdiff_lag6' ...
    '+ Fut19_firstdiff_lag7 + Fut19_firstdiff_lag8 + Fut19_firstdiff_lag9 + Fut19_firstdiff_lag10' ...
    '+ Fut19_firstdiff_lag11 + Fut19_firstdiff_lag12 + Spot_firstdiff_lag1 + Spot_firstdiff_lag2' ...
    '+ Spot_firstdiff_lag3 + Spot_firstdiff_lag4 + Spot_firstdiff_lag5 + Spot_firstdiff_lag6' ...
    '+ Spot_firstdiff_lag7 + Spot_firstdiff_lag8 + Spot_firstdiff_lag9 + Spot_firstdiff_lag10' ...
    '+ Spot_firstdiff_lag11 + Spot_firstdiff_lag12 + condVol'])

Pred19 = forecast(estMdl, height(residuals), 'Y0', ECM1.Spot_firstdiff);
Pred19_Vol = sqrt(Pred19);


Predicted_SpotFD_GARCH1 = [ECM_GARCH_1.Coefficients.Estimate(1) + ECM_GARCH_1.Coefficients.Estimate(3)*data.Fut19_firstdiff_lag1(1461:end) ...
    + ECM_GARCH_1.Coefficients.Estimate(4)*data.Fut19_firstdiff_lag2(1461:end) + ECM_GARCH_1.Coefficients.Estimate(5)*data.Fut19_firstdiff_lag3(1461:end) ...
    + ECM_GARCH_1.Coefficients.Estimate(6)*data.Fut19_firstdiff_lag4(1461:end) + ECM_GARCH_1.Coefficients.Estimate(7)*data.Fut19_firstdiff_lag5(1461:end) ...
    + ECM_GARCH_1.Coefficients.Estimate(8)*data.Fut19_firstdiff_lag6(1461:end) + ECM_GARCH_1.Coefficients.Estimate(9)*data.Fut19_firstdiff_lag7(1461:end) ...
    + ECM_GARCH_1.Coefficients.Estimate(10)*data.Fut19_firstdiff_lag8(1461:end) + ECM_GARCH_1.Coefficients.Estimate(11)*data.Fut19_firstdiff_lag9(1461:end) ...
    + ECM_GARCH_1.Coefficients.Estimate(12)*data.Fut19_firstdiff_lag10(1461:end) + ECM_GARCH_1.Coefficients.Estimate(13)*data.Fut19_firstdiff_lag11(1461:end) ...
    + ECM_GARCH_1.Coefficients.Estimate(14)*data.Fut19_firstdiff_lag12(1461:end) + ECM_GARCH_1.Coefficients.Estimate(15)*data.Spot_firstdiff_lag1(1461:end) ...
    + ECM_GARCH_1.Coefficients.Estimate(16)*data.Spot_firstdiff_lag2(1461:end) + ECM_GARCH_1.Coefficients.Estimate(17)*data.Spot_firstdiff_lag3(1461:end) ...
    + ECM_GARCH_1.Coefficients.Estimate(18)*data.Spot_firstdiff_lag4(1461:end) + ECM_GARCH_1.Coefficients.Estimate(19)*data.Spot_firstdiff_lag5(1461:end) ...
    + ECM_GARCH_1.Coefficients.Estimate(20)*data.Spot_firstdiff_lag6(1461:end) + ECM_GARCH_1.Coefficients.Estimate(21)*data.Spot_firstdiff_lag7(1461:end) ...
    + ECM_GARCH_1.Coefficients.Estimate(22)*data.Spot_firstdiff_lag8(1461:end) + ECM_GARCH_1.Coefficients.Estimate(23)*data.Spot_firstdiff_lag9(1461:end) ...
    + ECM_GARCH_1.Coefficients.Estimate(24)*data.Spot_firstdiff_lag10(1461:end) + ECM_GARCH_1.Coefficients.Estimate(25)*data.Spot_firstdiff_lag11(1461:end) ...
    + ECM_GARCH_1.Coefficients.Estimate(26)*data.Spot_firstdiff_lag12(1461:end) + ECM_GARCH_1.Coefficients.Estimate(2)*residuals.LaggedResiduals ...
    + ECM_GARCH_1.Coefficients.Estimate(27)*Pred19_Vol];

Predictions_GARCH1 = table();
Predictions_GARCH1.Date = data.Date(1460:end);
Predictions_GARCH1.Spot = data.Spot(1460:end);
Predictions_GARCH1.PredictedSpotFD = zeros(250,1);
Predictions_GARCH1.PredictedSpotFD(2:end) = Predicted_SpotFD_GARCH1;
Predictions_GARCH1.PredictedSpot = lagmatrix(Predictions_GARCH1.Spot, 1) + Predictions_GARCH1.PredictedSpotFD;
Predictions_GARCH1(1,:) = [];
Predictions_GARCH1.Spot = exp(Predictions_GARCH1.Spot);
Predictions_GARCH1.PredictedSpot = exp(Predictions_GARCH1.PredictedSpot);
Predictions_GARCH1(241:end,:) = [];

MSE_GARCH1 = (sum(Predictions_GARCH1.PredictedSpot - Predictions_GARCH1.Spot)^2)/(240);
RMSE_GARCH1 = sqrt(MSE_GARCH1)

ChangePred_GARCH1 = Predictions_GARCH1.PredictedSpot - lagmatrix(Predictions_GARCH1.PredictedSpot, 1);
ChangeSpot_GARCH1 = Predictions_GARCH1.Spot - lagmatrix(Predictions_GARCH1.Spot, 1);

SignPred_GARCH1 = table();
for i = 1:height(Predictions_GARCH1)
    if ChangePred_GARCH1(i) >= 0 & ChangeSpot_GARCH1(i) >= 0 
        SignPred_GARCH1.True(i) = 1;
    else
        SignPred_GARCH1.True(i) = 0;
    end
end

for i = 1:height(Predictions_GARCH1)
    if ChangePred_GARCH1(i) <= 0 & ChangeSpot_GARCH1(i) <= 0 
        SignPred_GARCH1.True1(i) = 1;
    else
        SignPred_GARCH1.True1(i) = 0;
    end
end
PercentCorrect_GARCH1 = (sum(SignPred_GARCH1.True) + sum(SignPred_GARCH1.True1)) / (height(Predictions_GARCH1));


% Finding residuals between Spot and lagged Fut18:
reg2 = table();
reg2.Date = data.Date(1:1214);
reg2.Spot = data.Spot(1:1214);
reg2.LagFut18 = lagmatrix(data.Fut18(1:1214),1);
reg2(1,:) = [];
model2 = fitlm(reg2, 'Spot ~ LagFut18')
reg2.residuals2 = model2.Residuals.Raw;

% Defining ECM variables - Spot and Fut18:
ECM2 = table();
ECM2.Date = reg2.Date;
ECM2.Spot_firstdiff = data.Spot_firstdiff(2:1214);
ECM2.LaggedResiduals = lagmatrix(reg2.residuals2,1);
% Lag variables - Fut:
ECM2.Fut18_firstdiff_lag1 = reg2.LagFut18 - lagmatrix(reg2.LagFut18,1);
ECM2.Fut18_firstdiff_lag2 = lagmatrix(ECM2.Fut18_firstdiff_lag1,1);
ECM2.Fut18_firstdiff_lag3 = lagmatrix(ECM2.Fut18_firstdiff_lag1,2);
ECM2.Fut18_firstdiff_lag4 = lagmatrix(ECM2.Fut18_firstdiff_lag1,3);
ECM2.Fut18_firstdiff_lag5 = lagmatrix(ECM2.Fut18_firstdiff_lag1,4);
ECM2.Fut18_firstdiff_lag6 = lagmatrix(ECM2.Fut18_firstdiff_lag1,5);
ECM2.Fut18_firstdiff_lag7 = lagmatrix(ECM2.Fut18_firstdiff_lag1,6);
ECM2.Fut18_firstdiff_lag8 = lagmatrix(ECM2.Fut18_firstdiff_lag1,7);
% Lag variables spot:
ECM2.Spot_firstdiff_lag1 = lagmatrix(ECM2.Spot_firstdiff,1);
ECM2.Spot_firstdiff_lag2 = lagmatrix(ECM2.Spot_firstdiff,2);
ECM2.Spot_firstdiff_lag3 = lagmatrix(ECM2.Spot_firstdiff,3);
ECM2.Spot_firstdiff_lag4 = lagmatrix(ECM2.Spot_firstdiff,4);
ECM2.Spot_firstdiff_lag5 = lagmatrix(ECM2.Spot_firstdiff,5);
ECM2.Spot_firstdiff_lag6 = lagmatrix(ECM2.Spot_firstdiff,6);
ECM2.Spot_firstdiff_lag7 = lagmatrix(ECM2.Spot_firstdiff,7);
ECM2.Spot_firstdiff_lag8 = lagmatrix(ECM2.Spot_firstdiff,8);

ECM_2 = fitlm(ECM2, ['Spot_firstdiff ~ LaggedResiduals + Fut18_firstdiff_lag1 + Fut18_firstdiff_lag2' ... 
    '+ Fut18_firstdiff_lag3 + Fut18_firstdiff_lag4 + Fut18_firstdiff_lag5 + Fut18_firstdiff_lag6' ...
    '+ Fut18_firstdiff_lag7 + Fut18_firstdiff_lag8 + Spot_firstdiff_lag1 + Spot_firstdiff_lag2' ...
    '+ Spot_firstdiff_lag3 + Spot_firstdiff_lag4 + Spot_firstdiff_lag5 + Spot_firstdiff_lag6' ...
    '+ Spot_firstdiff_lag7 + Spot_firstdiff_lag8'])

residuals_2 = table();
residuals_2.Date = data.Date;
residuals_2.LaggedResiduals = zeros(1709,1);

for i = 1214:1460
   reg = table();
   reg.Date = data.Date(1:i);
   reg.Spot = data.Spot(1:i);
   reg.LagFut18 = lagmatrix(data.Fut18(1:i),1);
   reg(1,:) = [];
   model = fitlm(reg, 'Spot ~ LagFut18');
   reg.residuals = model.Residuals.Raw; 
   residuals_2.LaggedResiduals(i+1) = reg.residuals(end);
end

residuals_2(1:1214,:) =[];
residuals_2(247:end,:) = [];
    
    
Predicted_SpotFD_ECM2 = [ECM_2.Coefficients.Estimate(1) + ECM_2.Coefficients.Estimate(3)*data.Fut18_firstdiff_lag1(1215:1460) ...
    + ECM_2.Coefficients.Estimate(4)*data.Fut18_firstdiff_lag2(1215:1460) + ECM_2.Coefficients.Estimate(5)*data.Fut18_firstdiff_lag3(1215:1460) ...
    + ECM_2.Coefficients.Estimate(6)*data.Fut18_firstdiff_lag4(1215:1460) + ECM_2.Coefficients.Estimate(7)*data.Fut18_firstdiff_lag5(1215:1460) ...
    + ECM_2.Coefficients.Estimate(8)*data.Fut18_firstdiff_lag6(1215:1460) + ECM_2.Coefficients.Estimate(9)*data.Fut18_firstdiff_lag7(1215:1460) ...
    + ECM_2.Coefficients.Estimate(10)*data.Fut18_firstdiff_lag8(1215:1460) + ECM_2.Coefficients.Estimate(11)*data.Spot_firstdiff_lag1(1215:1460) ...
    + ECM_2.Coefficients.Estimate(12)*data.Spot_firstdiff_lag2(1215:1460) + ECM_2.Coefficients.Estimate(13)*data.Spot_firstdiff_lag3(1215:1460) ...
    + ECM_2.Coefficients.Estimate(14)*data.Spot_firstdiff_lag4(1215:1460) + ECM_2.Coefficients.Estimate(15)*data.Spot_firstdiff_lag5(1215:1460) ...
    + ECM_2.Coefficients.Estimate(16)*data.Spot_firstdiff_lag6(1215:1460) + ECM_2.Coefficients.Estimate(17)*data.Spot_firstdiff_lag7(1215:1460) ...
    + ECM_2.Coefficients.Estimate(18)*data.Spot_firstdiff_lag8(1215:1460) + ECM_2.Coefficients.Estimate(2)*residuals_2.LaggedResiduals];

Predictions_ECM2 = table();
Predictions_ECM2.Date = data.Date(1214:1460);
Predictions_ECM2.Spot = data.Spot(1214:1460);
Predictions_ECM2.PredictedSpotFD = zeros(247,1);
Predictions_ECM2.PredictedSpotFD(2:end) = Predicted_SpotFD_ECM2;
Predictions_ECM2.PredictedSpot = lagmatrix(Predictions_ECM2.Spot, 1) + Predictions_ECM2.PredictedSpotFD;
Predictions_ECM2(1,:) = [];
Predictions_ECM2.Spot = exp(Predictions_ECM2.Spot);
Predictions_ECM2.PredictedSpot = exp(Predictions_ECM2.PredictedSpot);
Predictions_ECM2(242:end,:) = [];

MSE_ECM2 = (sum(Predictions_ECM2.PredictedSpot - Predictions_ECM2.Spot)^2)/(239);
RMSE_ECM2 = sqrt(MSE_ECM2)

ChangePred_ECM2 = Predictions_ECM2.PredictedSpot - lagmatrix(Predictions_ECM2.PredictedSpot, 1);
ChangeSpot_ECM2 = Predictions_ECM2.Spot - lagmatrix(Predictions_ECM2.Spot, 1);

SignPred_ECM2 = table();
for i = 1:height(Predictions_ECM2)
    if ChangePred_ECM2(i) >= 0 & ChangeSpot_ECM2(i) >= 0 
        SignPred_ECM2.True(i) = 1;
    else
        SignPred_ECM2.True(i) = 0;
    end
end

for i = 1:height(Predictions_ECM2)
    if ChangePred_ECM2(i) <= 0 & ChangeSpot_ECM2(i) <= 0 
        SignPred_ECM2.True1(i) = 1;
    else
        SignPred_ECM2.True1(i) = 0;
    end
end

PercentCorrect_ECM2 = (sum(SignPred_ECM2.True) + sum(SignPred_ECM2.True1)) / (height(Predictions_ECM2));


% ECM-GARCH: Spot and Fut18
GARCH2 = garch('GARCHlags', 1, 'ARCHlags', 1)
[estMdl, estParamCov, logL] = estimate(GARCH2, ECM2.Spot_firstdiff)
condVar_2 = infer(estMdl, ECM2.Spot_firstdiff);
ECM2.condVol = sqrt(condVar_2);

ECM_GARCH_2 = fitlm(ECM2, ['Spot_firstdiff ~ LaggedResiduals + Fut18_firstdiff_lag1 + Fut18_firstdiff_lag2' ... 
    '+ Fut18_firstdiff_lag3 + Fut18_firstdiff_lag4 + Fut18_firstdiff_lag5 + Fut18_firstdiff_lag6' ...
    '+ Fut18_firstdiff_lag7 + Fut18_firstdiff_lag8 + Spot_firstdiff_lag1 + Spot_firstdiff_lag2' ...
    '+ Spot_firstdiff_lag3 + Spot_firstdiff_lag4 + Spot_firstdiff_lag5 + Spot_firstdiff_lag6' ...
    '+ Spot_firstdiff_lag7 + Spot_firstdiff_lag8 + condVol'])

Pred18 = forecast(estMdl, height(residuals_2), 'Y0', ECM2.Spot_firstdiff);
Pred18_Vol = sqrt(Pred18);


Predicted_SpotFD_GARCH2 = [ECM_GARCH_2.Coefficients.Estimate(1) + ECM_GARCH_2.Coefficients.Estimate(3)*data.Fut18_firstdiff_lag1(1215:1460) ...
    + ECM_GARCH_2.Coefficients.Estimate(4)*data.Fut18_firstdiff_lag2(1215:1460) + ECM_GARCH_2.Coefficients.Estimate(5)*data.Fut18_firstdiff_lag3(1215:1460) ...
    + ECM_GARCH_2.Coefficients.Estimate(6)*data.Fut18_firstdiff_lag4(1215:1460) + ECM_GARCH_2.Coefficients.Estimate(7)*data.Fut18_firstdiff_lag5(1215:1460) ...
    + ECM_GARCH_2.Coefficients.Estimate(8)*data.Fut18_firstdiff_lag6(1215:1460) + ECM_GARCH_2.Coefficients.Estimate(9)*data.Fut18_firstdiff_lag7(1215:1460) ...
    + ECM_GARCH_2.Coefficients.Estimate(10)*data.Fut18_firstdiff_lag8(1215:1460) + ECM_GARCH_2.Coefficients.Estimate(11)*data.Spot_firstdiff_lag1(1215:1460) ...
    + ECM_GARCH_2.Coefficients.Estimate(12)*data.Spot_firstdiff_lag2(1215:1460) + ECM_GARCH_2.Coefficients.Estimate(13)*data.Spot_firstdiff_lag3(1215:1460) ...
    + ECM_GARCH_2.Coefficients.Estimate(14)*data.Spot_firstdiff_lag4(1215:1460) + ECM_GARCH_2.Coefficients.Estimate(15)*data.Spot_firstdiff_lag5(1215:1460) ...
    + ECM_GARCH_2.Coefficients.Estimate(16)*data.Spot_firstdiff_lag6(1215:1460) + ECM_GARCH_2.Coefficients.Estimate(17)*data.Spot_firstdiff_lag7(1215:1460) ...
    + ECM_GARCH_2.Coefficients.Estimate(18)*data.Spot_firstdiff_lag8(1215:1460) + ECM_GARCH_2.Coefficients.Estimate(2)*residuals_2.LaggedResiduals ... 
    + ECM_GARCH_2.Coefficients.Estimate(19)*Pred18_Vol];

Predictions_GARCH2 = table();
Predictions_GARCH2.Date = data.Date(1214:1460);
Predictions_GARCH2.Spot = data.Spot(1214:1460);
Predictions_GARCH2.PredictedSpotFD = zeros(247,1);
Predictions_GARCH2.PredictedSpotFD(2:end) = Predicted_SpotFD_GARCH2;
Predictions_GARCH2.PredictedSpot = lagmatrix(Predictions_GARCH2.Spot, 1) + Predictions_GARCH2.PredictedSpotFD;
Predictions_GARCH2(1,:) = [];
Predictions_GARCH2.Spot = exp(Predictions_GARCH2.Spot);
Predictions_GARCH2.PredictedSpot = exp(Predictions_GARCH2.PredictedSpot);
Predictions_GARCH2(242:end,:) = [];

MSE_GARCH2 = (sum(Predictions_GARCH2.PredictedSpot - Predictions_GARCH2.Spot)^2)/(239);
RMSE_GARCH2 = sqrt(MSE_GARCH2)

ChangePred_GARCH2 = Predictions_GARCH2.PredictedSpot - lagmatrix(Predictions_GARCH2.PredictedSpot, 1);
ChangeSpot_GARCH2 = Predictions_GARCH2.Spot - lagmatrix(Predictions_GARCH2.Spot, 1);

SignPred_GARCH2 = table();
for i = 1:height(Predictions_GARCH2)
    if ChangePred_GARCH2(i) >= 0 & ChangeSpot_GARCH2(i) >= 0 
        SignPred_GARCH2.True(i) = 1;
    else
        SignPred_GARCH2.True(i) = 0;
    end
end

for i = 1:239
    if ChangePred_GARCH2(i) <= 0 & ChangeSpot_GARCH2(i) <= 0 
        SignPred_GARCH2.True1(i) = 1;
    else
        SignPred_GARCH2.True1(i) = 0;
    end
end

PercentCorrect_GARCH2 = (sum(SignPred_GARCH2.True) + sum(SignPred_GARCH2.True1)) / (height(Predictions_GARCH2));


% Finding residuals between Spot and lagged Fut17:
reg3 = table();
reg3.Date = data.Date(1:990);
reg3.Spot = data.Spot(1:990);
reg3.LagFut17 = lagmatrix(data.Fut17(1:990),1);
reg3(1,:) = [];
model3 = fitlm(reg3, 'Spot ~ LagFut17')
reg3.residuals3 = model3.Residuals.Raw;

% Defining ECM variables - Spot and Fut18:
ECM3 = table();
ECM3.Date = reg3.Date;
ECM3.Spot_firstdiff = data.Spot_firstdiff(2:990);
ECM3.LaggedResiduals = lagmatrix(reg3.residuals3,1);
%Lag variables Fut17:
ECM3.Fut17_firstdiff_lag1 = reg3.LagFut17 - lagmatrix(reg3.LagFut17,1);
ECM3.Fut17_firstdiff_lag2 = lagmatrix(ECM3.Fut17_firstdiff_lag1,1);
ECM3.Fut17_firstdiff_lag3 = lagmatrix(ECM3.Fut17_firstdiff_lag1,2);
ECM3.Fut17_firstdiff_lag4 = lagmatrix(ECM3.Fut17_firstdiff_lag1,3);
ECM3.Fut17_firstdiff_lag5 = lagmatrix(ECM3.Fut17_firstdiff_lag1,4);
ECM3.Fut17_firstdiff_lag6 = lagmatrix(ECM3.Fut17_firstdiff_lag1,5);
ECM3.Fut17_firstdiff_lag7 = lagmatrix(ECM3.Fut17_firstdiff_lag1,6);
ECM3.Fut17_firstdiff_lag8 = lagmatrix(ECM3.Fut17_firstdiff_lag1,7);
ECM3.Fut17_firstdiff_lag9 = lagmatrix(ECM3.Fut17_firstdiff_lag1,8);
ECM3.Fut17_firstdiff_lag10 = lagmatrix(ECM3.Fut17_firstdiff_lag1,9);
ECM3.Fut17_firstdiff_lag11 = lagmatrix(ECM3.Fut17_firstdiff_lag1,10);
ECM3.Fut17_firstdiff_lag12 = lagmatrix(ECM3.Fut17_firstdiff_lag1,11);

%Lag variables spot:
ECM3.Spot_firstdiff_lag1 = lagmatrix(ECM3.Spot_firstdiff,1);
ECM3.Spot_firstdiff_lag2 = lagmatrix(ECM3.Spot_firstdiff,2);
ECM3.Spot_firstdiff_lag3 = lagmatrix(ECM3.Spot_firstdiff,3);
ECM3.Spot_firstdiff_lag4 = lagmatrix(ECM3.Spot_firstdiff,4);
ECM3.Spot_firstdiff_lag5 = lagmatrix(ECM3.Spot_firstdiff,5);
ECM3.Spot_firstdiff_lag6 = lagmatrix(ECM3.Spot_firstdiff,6);
ECM3.Spot_firstdiff_lag7 = lagmatrix(ECM3.Spot_firstdiff,7);
ECM3.Spot_firstdiff_lag8 = lagmatrix(ECM3.Spot_firstdiff,8);
ECM3.Spot_firstdiff_lag9 = lagmatrix(ECM3.Spot_firstdiff,9);
ECM3.Spot_firstdiff_lag10 = lagmatrix(ECM3.Spot_firstdiff,10);
ECM3.Spot_firstdiff_lag11 = lagmatrix(ECM3.Spot_firstdiff,11);
ECM3.Spot_firstdiff_lag12 = lagmatrix(ECM3.Spot_firstdiff,12);

ECM_3 = fitlm(ECM3, ['Spot_firstdiff ~ LaggedResiduals + Fut17_firstdiff_lag1 + Fut17_firstdiff_lag2' ... 
    '+ Fut17_firstdiff_lag3 + Fut17_firstdiff_lag4 + Fut17_firstdiff_lag5 + Fut17_firstdiff_lag6' ...
    '+ Fut17_firstdiff_lag7 + Fut17_firstdiff_lag8 + Fut17_firstdiff_lag9 + Fut17_firstdiff_lag10' ...
    '+ Fut17_firstdiff_lag11 + Fut17_firstdiff_lag12 + Spot_firstdiff_lag1 + Spot_firstdiff_lag2' ...
    '+ Spot_firstdiff_lag3 + Spot_firstdiff_lag4 + Spot_firstdiff_lag5 + Spot_firstdiff_lag6' ...
    '+ Spot_firstdiff_lag7 + Spot_firstdiff_lag8 + Spot_firstdiff_lag9 + Spot_firstdiff_lag10' ...
    '+ Spot_firstdiff_lag11 + Spot_firstdiff_lag12'])

residuals_3 = table();
residuals_3.Date = data.Date;
residuals_3.LaggedResiduals = zeros(1709,1);

for i = 990:1214
   reg = table();
   reg.Date = data.Date(1:i);
   reg.Spot = data.Spot(1:i);
   reg.LagFut17 = lagmatrix(data.Fut17(1:i),1);
   reg(1,:) = [];
   model = fitlm(reg, 'Spot ~ LagFut17');
   reg.residuals = model.Residuals.Raw; 
   residuals_3.LaggedResiduals(i+1) = reg.residuals(end);
end

residuals_3(1:990,:) =[];
residuals_3(225:end,:) = [];
    
    
Predicted_SpotFD_ECM3 = [ECM_3.Coefficients.Estimate(1) + ECM_3.Coefficients.Estimate(3)*data.Fut17_firstdiff_lag1(991:1214) ...
    + ECM_3.Coefficients.Estimate(4)*data.Fut17_firstdiff_lag2(991:1214) + ECM_3.Coefficients.Estimate(5)*data.Fut17_firstdiff_lag3(991:1214) ...
    + ECM_3.Coefficients.Estimate(6)*data.Fut17_firstdiff_lag4(991:1214) + ECM_3.Coefficients.Estimate(7)*data.Fut17_firstdiff_lag5(991:1214) ...
    + ECM_3.Coefficients.Estimate(8)*data.Fut17_firstdiff_lag6(991:1214) + ECM_3.Coefficients.Estimate(9)*data.Fut17_firstdiff_lag7(991:1214) ...
    + ECM_3.Coefficients.Estimate(10)*data.Fut17_firstdiff_lag8(991:1214) + ECM_3.Coefficients.Estimate(11)*data.Fut17_firstdiff_lag9(991:1214) ...
    + ECM_3.Coefficients.Estimate(12)*data.Fut17_firstdiff_lag10(991:1214) + ECM_3.Coefficients.Estimate(13)*data.Fut17_firstdiff_lag11(991:1214) ...
    + ECM_3.Coefficients.Estimate(14)*data.Fut17_firstdiff_lag12(991:1214) + ECM_3.Coefficients.Estimate(15)*data.Spot_firstdiff_lag1(991:1214) ...
    + ECM_3.Coefficients.Estimate(16)*data.Spot_firstdiff_lag2(991:1214) + ECM_3.Coefficients.Estimate(17)*data.Spot_firstdiff_lag3(991:1214) ...
    + ECM_3.Coefficients.Estimate(18)*data.Spot_firstdiff_lag4(991:1214) + ECM_3.Coefficients.Estimate(19)*data.Spot_firstdiff_lag5(991:1214) ...
    + ECM_3.Coefficients.Estimate(20)*data.Spot_firstdiff_lag6(991:1214) + ECM_3.Coefficients.Estimate(21)*data.Spot_firstdiff_lag7(991:1214) ...
    + ECM_3.Coefficients.Estimate(22)*data.Spot_firstdiff_lag8(991:1214) + ECM_3.Coefficients.Estimate(23)*data.Spot_firstdiff_lag9(991:1214) ...
    + ECM_3.Coefficients.Estimate(24)*data.Spot_firstdiff_lag10(991:1214) + ECM_3.Coefficients.Estimate(25)*data.Spot_firstdiff_lag11(991:1214) ...
    + ECM_3.Coefficients.Estimate(26)*data.Spot_firstdiff_lag12(991:1214) + ECM_3.Coefficients.Estimate(2)*residuals_3.LaggedResiduals];
   

Predictions_ECM3 = table();
Predictions_ECM3.Date = data.Date(990:1214);
Predictions_ECM3.Spot = data.Spot(990:1214);
Predictions_ECM3.PredictedSpotFD = zeros(225,1);
Predictions_ECM3.PredictedSpotFD(2:end) = Predicted_SpotFD_ECM3;
Predictions_ECM3.PredictedSpot = lagmatrix(Predictions_ECM3.Spot, 1) + Predictions_ECM3.PredictedSpotFD;
Predictions_ECM3(1,:) = [];
Predictions_ECM3.Spot = exp(Predictions_ECM3.Spot);
Predictions_ECM3.PredictedSpot = exp(Predictions_ECM3.PredictedSpot);
Predictions_ECM3(219:end,:) = [];

MSE_ECM3 = (sum(Predictions_ECM3.PredictedSpot - Predictions_ECM3.Spot)^2)/(218);
RMSE_ECM3 = sqrt(MSE_ECM3)

ChangePred_ECM3 = Predictions_ECM3.PredictedSpot - lagmatrix(Predictions_ECM3.PredictedSpot, 1);
ChangeSpot_ECM3 = Predictions_ECM3.Spot - lagmatrix(Predictions_ECM3.Spot, 1);

SignPred_ECM3 = table();
for i = 1:height(Predictions_ECM3)
    if ChangePred_ECM3(i) >= 0 & ChangeSpot_ECM3(i) >= 0 
        SignPred_ECM3.True(i) = 1;
    else
        SignPred_ECM3.True(i) = 0;
    end
end

for i = 1:height(Predictions_ECM3)
    if ChangePred_ECM3(i) <= 0 & ChangeSpot_ECM3(i) <= 0 
        SignPred_ECM3.True1(i) = 1;
    else
        SignPred_ECM3.True1(i) = 0;
    end
end

PercentCorrect_ECM3 = (sum(SignPred_ECM3.True) + sum(SignPred_ECM3.True1)) / (height(Predictions_ECM3));


% ECM-GARCH: Spot and Fut17
GARCH3 = garch('GARCHlags', 1, 'ARCHlags', 1)
[estMdl, estParamCov, logL] = estimate(GARCH3, ECM3.Spot_firstdiff)
condVar_3 = infer(estMdl, ECM3.Spot_firstdiff);
ECM3.condVol = sqrt(condVar_3);

ECM_GARCH_3 = fitlm(ECM3, ['Spot_firstdiff ~ LaggedResiduals + Fut17_firstdiff_lag1 + Fut17_firstdiff_lag2' ... 
    '+ Fut17_firstdiff_lag3 + Fut17_firstdiff_lag4 + Fut17_firstdiff_lag5 + Fut17_firstdiff_lag6' ...
    '+ Fut17_firstdiff_lag7 + Fut17_firstdiff_lag8 + Fut17_firstdiff_lag9 + Fut17_firstdiff_lag10' ...
    '+ Fut17_firstdiff_lag11 + Fut17_firstdiff_lag12 + Spot_firstdiff_lag1 + Spot_firstdiff_lag2' ...
    '+ Spot_firstdiff_lag3 + Spot_firstdiff_lag4 + Spot_firstdiff_lag5 + Spot_firstdiff_lag6' ...
    '+ Spot_firstdiff_lag7 + Spot_firstdiff_lag8 + Spot_firstdiff_lag9 + Spot_firstdiff_lag10' ...
    '+ Spot_firstdiff_lag11 + Spot_firstdiff_lag12 + condVol'])

Pred17 = forecast(estMdl, height(residuals_3), 'Y0', ECM3.Spot_firstdiff);
Pred17_Vol = sqrt(Pred17);


Predicted_SpotFD_GARCH3 = [ECM_GARCH_3.Coefficients.Estimate(1) + ECM_GARCH_3.Coefficients.Estimate(3)*data.Fut17_firstdiff_lag1(991:1214) ...
    + ECM_GARCH_3.Coefficients.Estimate(4)*data.Fut17_firstdiff_lag2(991:1214) + ECM_GARCH_3.Coefficients.Estimate(5)*data.Fut17_firstdiff_lag3(991:1214) ...
    + ECM_GARCH_3.Coefficients.Estimate(6)*data.Fut17_firstdiff_lag4(991:1214) + ECM_GARCH_3.Coefficients.Estimate(7)*data.Fut17_firstdiff_lag5(991:1214) ...
    + ECM_GARCH_3.Coefficients.Estimate(8)*data.Fut17_firstdiff_lag6(991:1214) + ECM_GARCH_3.Coefficients.Estimate(9)*data.Fut17_firstdiff_lag7(991:1214) ...
    + ECM_GARCH_3.Coefficients.Estimate(10)*data.Fut17_firstdiff_lag8(991:1214) + ECM_GARCH_3.Coefficients.Estimate(11)*data.Fut17_firstdiff_lag9(991:1214) ...
    + ECM_GARCH_3.Coefficients.Estimate(12)*data.Fut17_firstdiff_lag10(991:1214) + ECM_GARCH_3.Coefficients.Estimate(13)*data.Fut17_firstdiff_lag11(991:1214) ...
    + ECM_GARCH_3.Coefficients.Estimate(14)*data.Fut17_firstdiff_lag12(991:1214) + ECM_GARCH_3.Coefficients.Estimate(15)*data.Spot_firstdiff_lag1(991:1214) ...
    + ECM_GARCH_3.Coefficients.Estimate(16)*data.Spot_firstdiff_lag2(991:1214) + ECM_GARCH_3.Coefficients.Estimate(17)*data.Spot_firstdiff_lag3(991:1214) ...
    + ECM_GARCH_3.Coefficients.Estimate(18)*data.Spot_firstdiff_lag4(991:1214) + ECM_GARCH_3.Coefficients.Estimate(19)*data.Spot_firstdiff_lag5(991:1214) ...
    + ECM_GARCH_3.Coefficients.Estimate(20)*data.Spot_firstdiff_lag6(991:1214) + ECM_GARCH_3.Coefficients.Estimate(21)*data.Spot_firstdiff_lag7(991:1214) ...
    + ECM_GARCH_3.Coefficients.Estimate(22)*data.Spot_firstdiff_lag8(991:1214) + ECM_GARCH_3.Coefficients.Estimate(23)*data.Spot_firstdiff_lag9(991:1214) ...
    + ECM_GARCH_3.Coefficients.Estimate(24)*data.Spot_firstdiff_lag10(991:1214) + ECM_GARCH_3.Coefficients.Estimate(25)*data.Spot_firstdiff_lag11(991:1214) ...
    + ECM_GARCH_3.Coefficients.Estimate(26)*data.Spot_firstdiff_lag12(991:1214) + ECM_GARCH_3.Coefficients.Estimate(2)*residuals_3.LaggedResiduals ...
    + ECM_GARCH_3.Coefficients.Estimate(27)*Pred17_Vol];

Predictions_GARCH3 = table();
Predictions_GARCH3.Date = data.Date(990:1214);
Predictions_GARCH3.Spot = data.Spot(990:1214);
Predictions_GARCH3.PredictedSpotFD = zeros(225,1);
Predictions_GARCH3.PredictedSpotFD(2:end) = Predicted_SpotFD_GARCH3;
Predictions_GARCH3.PredictedSpot = lagmatrix(Predictions_GARCH3.Spot, 1) + Predictions_GARCH3.PredictedSpotFD;
Predictions_GARCH3(1,:) = [];
Predictions_GARCH3.Spot = exp(Predictions_GARCH3.Spot);
Predictions_GARCH3.PredictedSpot = exp(Predictions_GARCH3.PredictedSpot);
Predictions_GARCH3(219:end,:) = [];

MSE_GARCH3 = (sum(Predictions_GARCH3.PredictedSpot - Predictions_GARCH3.Spot)^2)/(218);
RMSE_GARCH3 = sqrt(MSE_GARCH3)

ChangePred_GARCH3 = Predictions_GARCH3.PredictedSpot - lagmatrix(Predictions_GARCH3.PredictedSpot, 1);
ChangeSpot_GARCH3 = Predictions_GARCH3.Spot - lagmatrix(Predictions_GARCH3.Spot, 1);

SignPred_GARCH3 = table();
for i = 1:height(Predictions_GARCH3)
    if ChangePred_GARCH3(i) >= 0 & ChangeSpot_GARCH3(i) >= 0 
        SignPred_GARCH3.True(i) = 1;
    else
        SignPred_GARCH3.True(i) = 0;
    end
end

for i = 1:height(Predictions_GARCH3)
    if ChangePred_GARCH3(i) <= 0 & ChangeSpot_GARCH3(i) <= 0 
        SignPred_GARCH3.True1(i) = 1;
    else
        SignPred_GARCH3.True1(i) = 0;
    end
end
PercentCorrect_GARCH3 = (sum(SignPred_GARCH3.True) + sum(SignPred_GARCH3.True1)) / (height(Predictions_GARCH3));


% Finding residuals between Spot and lagged Fut16:
reg4 = table();
reg4.Date = data.Date(1:739);
reg4.Spot = data.Spot(1:739);
reg4.LagFut16 = lagmatrix(data.Fut16(1:739),1);
reg4(1,:) = [];
model4 = fitlm(reg4, 'Spot ~ LagFut16')
reg4.residuals4 = model4.Residuals.Raw;

% Defining ECM variables - Spot and Fut16:
ECM4 = table();
ECM4.Date = reg4.Date;
ECM4.Spot_firstdiff = data.Spot_firstdiff(2:739);
ECM4.LaggedResiduals = lagmatrix(reg4.residuals4,1);
%Lag variables - Fut:
ECM4.Fut16_firstdiff_lag1 = reg4.LagFut16 - lagmatrix(reg4.LagFut16,1);
ECM4.Fut16_firstdiff_lag2 = lagmatrix(ECM4.Fut16_firstdiff_lag1,1);
ECM4.Fut16_firstdiff_lag3 = lagmatrix(ECM4.Fut16_firstdiff_lag1,2);
ECM4.Fut16_firstdiff_lag4 = lagmatrix(ECM4.Fut16_firstdiff_lag1,3);
ECM4.Fut16_firstdiff_lag5 = lagmatrix(ECM4.Fut16_firstdiff_lag1,4);
ECM4.Fut16_firstdiff_lag6 = lagmatrix(ECM4.Fut16_firstdiff_lag1,5);
ECM4.Fut16_firstdiff_lag7 = lagmatrix(ECM4.Fut16_firstdiff_lag1,6);


%Lag variables spot:
ECM4.Spot_firstdiff_lag1 = lagmatrix(ECM4.Spot_firstdiff,1);
ECM4.Spot_firstdiff_lag2 = lagmatrix(ECM4.Spot_firstdiff,2);
ECM4.Spot_firstdiff_lag3 = lagmatrix(ECM4.Spot_firstdiff,3);
ECM4.Spot_firstdiff_lag4 = lagmatrix(ECM4.Spot_firstdiff,4);
ECM4.Spot_firstdiff_lag5 = lagmatrix(ECM4.Spot_firstdiff,5);
ECM4.Spot_firstdiff_lag6 = lagmatrix(ECM4.Spot_firstdiff,6);
ECM4.Spot_firstdiff_lag7 = lagmatrix(ECM4.Spot_firstdiff,7);



ECM_4 = fitlm(ECM4, ['Spot_firstdiff ~ LaggedResiduals + Fut16_firstdiff_lag1 + Fut16_firstdiff_lag2' ... 
    '+ Fut16_firstdiff_lag3 + Fut16_firstdiff_lag4 + Fut16_firstdiff_lag5 + Fut16_firstdiff_lag6' ...
    '+ Fut16_firstdiff_lag7 + Spot_firstdiff_lag1 + Spot_firstdiff_lag2' ...
    '+ Spot_firstdiff_lag3 + Spot_firstdiff_lag4 + Spot_firstdiff_lag5 + Spot_firstdiff_lag6' ...
    '+ Spot_firstdiff_lag7'])

residuals_4 = table();
residuals_4.Date = data.Date;
residuals_4.LaggedResiduals = zeros(1709,1);

for i = 739:990
   reg = table();
   reg.Date = data.Date(1:i);
   reg.Spot = data.Spot(1:i);
   reg.LagFut16 = lagmatrix(data.Fut16(1:i),1);
   reg(1,:) = [];
   model4 = fitlm(reg, 'Spot ~ LagFut16');
   reg.residuals = model4.Residuals.Raw; 
   residuals_4.LaggedResiduals(i+1) = reg.residuals(end);
end

residuals_4(1:739,:) =[];
residuals_4(252:end,:) = [];

Predicted_SpotFD_ECM4 = [ECM_4.Coefficients.Estimate(1) + ECM_4.Coefficients.Estimate(3)*data.Fut16_firstdiff_lag1(740:990) ...
    + ECM_4.Coefficients.Estimate(4)*data.Fut16_firstdiff_lag2(740:990) + ECM_4.Coefficients.Estimate(5)*data.Fut16_firstdiff_lag3(740:990) ...
    + ECM_4.Coefficients.Estimate(6)*data.Fut16_firstdiff_lag4(740:990) + ECM_4.Coefficients.Estimate(7)*data.Fut16_firstdiff_lag5(740:990) ...
    + ECM_4.Coefficients.Estimate(8)*data.Fut16_firstdiff_lag6(740:990) + ECM_4.Coefficients.Estimate(9)*data.Fut16_firstdiff_lag7(740:990) ...
    + ECM_4.Coefficients.Estimate(10)*data.Spot_firstdiff_lag1(740:990) + ECM_4.Coefficients.Estimate(11)*data.Spot_firstdiff_lag2(740:990) ...
    + ECM_4.Coefficients.Estimate(12)*data.Spot_firstdiff_lag3(740:990) + ECM_4.Coefficients.Estimate(13)*data.Spot_firstdiff_lag4(740:990) ...
    + ECM_4.Coefficients.Estimate(14)*data.Spot_firstdiff_lag5(740:990) + ECM_4.Coefficients.Estimate(15)*data.Spot_firstdiff_lag6(740:990) ...
    + ECM_4.Coefficients.Estimate(16)*data.Spot_firstdiff_lag7(740:990) + ECM_4.Coefficients.Estimate(2)*residuals_4.LaggedResiduals];

Predictions_ECM4 = table();
Predictions_ECM4.Date = data.Date(739:990);
Predictions_ECM4.Spot = data.Spot(739:990);
Predictions_ECM4.PredictedSpotFD = zeros(252,1);
Predictions_ECM4.PredictedSpotFD(2:end) = Predicted_SpotFD_ECM4;
Predictions_ECM4.PredictedSpot = lagmatrix(Predictions_ECM4.Spot, 1) + Predictions_ECM4.PredictedSpotFD;
Predictions_ECM4(1,:) = [];
Predictions_ECM4.Spot = exp(Predictions_ECM4.Spot);
Predictions_ECM4.PredictedSpot = exp(Predictions_ECM4.PredictedSpot);
Predictions_ECM4(246:end,:) = [];

MSE_ECM4 = (sum(Predictions_ECM4.PredictedSpot - Predictions_ECM4.Spot)^2)/(245);
RMSE_ECM4 = sqrt(MSE_ECM4)

ChangePred_ECM4 = Predictions_ECM4.PredictedSpot - lagmatrix(Predictions_ECM4.PredictedSpot, 1);
ChangeSpot_ECM4 = Predictions_ECM4.Spot - lagmatrix(Predictions_ECM4.Spot, 1);

SignPred_ECM4 = table();
for i = 1:height(Predictions_ECM4)
    if ChangePred_ECM4(i) >= 0 & ChangeSpot_ECM4(i) >= 0 
        SignPred_ECM4.True(i) = 1;
    else
        SignPred_ECM4.True(i) = 0;
    end
end

for i = 1:height(Predictions_ECM4)
    if ChangePred_ECM4(i) <= 0 & ChangeSpot_ECM4(i) <= 0 
        SignPred_ECM4.True1(i) = 1;
    else
        SignPred_ECM4.True1(i) = 0;
    end
end

PercentCorrect_ECM4 = (sum(SignPred_ECM4.True) + sum(SignPred_ECM4.True1)) / (height(Predictions_ECM4));


% ECM-GARCH: Spot and Fut16
GARCH4 = garch('GARCHlags', 1, 'ARCHlags', 1)
[estMdl, estParamCov, logL] = estimate(GARCH4, ECM4.Spot_firstdiff)
condVar_4 = infer(estMdl, ECM4.Spot_firstdiff);
ECM4.condVol = sqrt(condVar_4);

ECM_GARCH_4 = fitlm(ECM4, ['Spot_firstdiff ~ LaggedResiduals + Fut16_firstdiff_lag1 + Fut16_firstdiff_lag2' ... 
    '+ Fut16_firstdiff_lag3 + Fut16_firstdiff_lag4 + Fut16_firstdiff_lag5 + Fut16_firstdiff_lag6' ...
    '+ Fut16_firstdiff_lag7 + Spot_firstdiff_lag1 + Spot_firstdiff_lag2' ...
    '+ Spot_firstdiff_lag3 + Spot_firstdiff_lag4 + Spot_firstdiff_lag5 + Spot_firstdiff_lag6' ...
    '+ Spot_firstdiff_lag7 + condVol'])

Pred16 = forecast(estMdl, height(residuals_4), 'Y0', ECM4.Spot_firstdiff);
Pred16_Vol = sqrt(Pred16);


Predicted_SpotFD_GARCH4 = [ECM_GARCH_4.Coefficients.Estimate(1) + ECM_GARCH_4.Coefficients.Estimate(3)*data.Fut16_firstdiff_lag1(740:990) ...
    + ECM_GARCH_4.Coefficients.Estimate(4)*data.Fut16_firstdiff_lag2(740:990) + ECM_GARCH_4.Coefficients.Estimate(5)*data.Fut16_firstdiff_lag3(740:990) ...
    + ECM_GARCH_4.Coefficients.Estimate(6)*data.Fut16_firstdiff_lag4(740:990) + ECM_GARCH_4.Coefficients.Estimate(7)*data.Fut16_firstdiff_lag5(740:990) ...
    + ECM_GARCH_4.Coefficients.Estimate(8)*data.Fut16_firstdiff_lag6(740:990) + ECM_GARCH_4.Coefficients.Estimate(9)*data.Fut16_firstdiff_lag7(740:990) ...
    + ECM_GARCH_4.Coefficients.Estimate(10)*data.Spot_firstdiff_lag1(740:990) + ECM_GARCH_4.Coefficients.Estimate(11)*data.Spot_firstdiff_lag2(740:990) ...
    + ECM_GARCH_4.Coefficients.Estimate(12)*data.Spot_firstdiff_lag3(740:990) + ECM_GARCH_4.Coefficients.Estimate(13)*data.Spot_firstdiff_lag4(740:990) ...
    + ECM_GARCH_4.Coefficients.Estimate(14)*data.Spot_firstdiff_lag5(740:990) + ECM_GARCH_4.Coefficients.Estimate(15)*data.Spot_firstdiff_lag6(740:990) ...
    + ECM_GARCH_4.Coefficients.Estimate(16)*data.Spot_firstdiff_lag7(740:990) + ECM_GARCH_4.Coefficients.Estimate(2)*residuals_4.LaggedResiduals ...
    + ECM_GARCH_4.Coefficients.Estimate(17)*Pred16_Vol];

Predictions_GARCH4 = table();
Predictions_GARCH4.Date = data.Date(739:990);
Predictions_GARCH4.Spot = data.Spot(739:990);
Predictions_GARCH4.PredictedSpotFD = zeros(252,1);
Predictions_GARCH4.PredictedSpotFD(2:end) = Predicted_SpotFD_GARCH4;
Predictions_GARCH4.PredictedSpot = lagmatrix(Predictions_GARCH4.Spot, 1) + Predictions_GARCH4.PredictedSpotFD;
Predictions_GARCH4(1,:) = [];
Predictions_GARCH4.Spot = exp(Predictions_GARCH4.Spot);
Predictions_GARCH4.PredictedSpot = exp(Predictions_GARCH4.PredictedSpot);
Predictions_GARCH4(246:end,:) = [];

MSE_GARCH4 = (sum(Predictions_GARCH4.PredictedSpot - Predictions_GARCH4.Spot)^2)/(245);
RMSE_GARCH4 = sqrt(MSE_GARCH4)

ChangePred_GARCH4 = Predictions_GARCH4.PredictedSpot - lagmatrix(Predictions_GARCH4.PredictedSpot, 1);
ChangeSpot_GARCH4 = Predictions_GARCH4.Spot - lagmatrix(Predictions_GARCH4.Spot, 1);

SignPred_GARCH4 = table();
for i = 1:height(Predictions_GARCH4)
    if ChangePred_GARCH4(i) >= 0 & ChangeSpot_GARCH4(i) >= 0 
        SignPred_GARCH4.True(i) = 1;
    else
        SignPred_GARCH4.True(i) = 0;
    end
end

for i = 1:height(Predictions_GARCH4)
    if ChangePred_GARCH4(i) <= 0 & ChangeSpot_GARCH4(i) <= 0 
        SignPred_GARCH4.True1(i) = 1;
    else
        SignPred_GARCH4.True1(i) = 0;
    end
end
PercentCorrect_GARCH4 = (sum(SignPred_GARCH4.True) + sum(SignPred_GARCH4.True1)) / (height(Predictions_GARCH4));

% Finding residuals between Spot and lagged Fut16:
reg5 = table();
reg5.Date = data.Date(1:493);
reg5.Spot = data.Spot(1:493);
reg5.LagFut15 = lagmatrix(data.Fut15(1:493),1);
reg5(1,:) = [];
model5 = fitlm(reg5, 'Spot ~ LagFut15')
reg5.residuals5 = model5.Residuals.Raw;

% Defining ECM variables - Spot and Fut15:
ECM5 = table();
ECM5.Date = reg5.Date;
ECM5.Spot_firstdiff = data.Spot_firstdiff(2:493);
ECM5.LaggedResiduals = lagmatrix(reg5.residuals5,1);
%Lag variables - Fut:
ECM5.Fut15_firstdiff_lag1 = reg5.LagFut15 - lagmatrix(reg5.LagFut15,1);
ECM5.Fut15_firstdiff_lag2 = lagmatrix(ECM5.Fut15_firstdiff_lag1,1);
ECM5.Fut15_firstdiff_lag3 = lagmatrix(ECM5.Fut15_firstdiff_lag1,2);
ECM5.Fut15_firstdiff_lag4 = lagmatrix(ECM5.Fut15_firstdiff_lag1,3);
ECM5.Fut15_firstdiff_lag5 = lagmatrix(ECM5.Fut15_firstdiff_lag1,4);
ECM5.Fut15_firstdiff_lag6 = lagmatrix(ECM5.Fut15_firstdiff_lag1,5);
ECM5.Fut15_firstdiff_lag7 = lagmatrix(ECM5.Fut15_firstdiff_lag1,6);
ECM5.Fut15_firstdiff_lag8 = lagmatrix(ECM5.Fut15_firstdiff_lag1,7);
ECM5.Fut15_firstdiff_lag9 = lagmatrix(ECM5.Fut15_firstdiff_lag1,8);
ECM5.Fut15_firstdiff_lag10 = lagmatrix(ECM5.Fut15_firstdiff_lag1,9);
%Lag variables spot:
ECM5.Spot_firstdiff_lag1 = lagmatrix(ECM5.Spot_firstdiff,1);
ECM5.Spot_firstdiff_lag2 = lagmatrix(ECM5.Spot_firstdiff,2);
ECM5.Spot_firstdiff_lag3 = lagmatrix(ECM5.Spot_firstdiff,3);
ECM5.Spot_firstdiff_lag4 = lagmatrix(ECM5.Spot_firstdiff,4);
ECM5.Spot_firstdiff_lag5 = lagmatrix(ECM5.Spot_firstdiff,5);
ECM5.Spot_firstdiff_lag6 = lagmatrix(ECM5.Spot_firstdiff,6);
ECM5.Spot_firstdiff_lag7 = lagmatrix(ECM5.Spot_firstdiff,7);
ECM5.Spot_firstdiff_lag8 = lagmatrix(ECM5.Spot_firstdiff,8);
ECM5.Spot_firstdiff_lag9 = lagmatrix(ECM5.Spot_firstdiff,9);
ECM5.Spot_firstdiff_lag10 = lagmatrix(ECM5.Spot_firstdiff,10);

ECM_5 = fitlm(ECM5, ['Spot_firstdiff ~ LaggedResiduals + Fut15_firstdiff_lag1 + Fut15_firstdiff_lag2' ... 
    '+ Fut15_firstdiff_lag3 + Fut15_firstdiff_lag4 + Fut15_firstdiff_lag5 + Fut15_firstdiff_lag6' ...
    '+ Fut15_firstdiff_lag7 + Fut15_firstdiff_lag8 + Fut15_firstdiff_lag9 + Fut15_firstdiff_lag10' ...
    '+ Spot_firstdiff_lag1 + Spot_firstdiff_lag2' ...
    '+ Spot_firstdiff_lag3 + Spot_firstdiff_lag4 + Spot_firstdiff_lag5 + Spot_firstdiff_lag6' ...
    '+ Spot_firstdiff_lag7 + Spot_firstdiff_lag8 + Spot_firstdiff_lag9 + Spot_firstdiff_lag10'])

residuals_5 = table();
residuals_5.Date = data.Date;
residuals_5.LaggedResiduals = zeros(1709,1);

for i = 493:739
   reg = table();
   reg.Date = data.Date(1:i);
   reg.Spot = data.Spot(1:i);
   reg.LagFut15 = lagmatrix(data.Fut15(1:i),1);
   reg(1,:) = [];
   model5 = fitlm(reg, 'Spot ~ LagFut15');
   reg.residuals = model5.Residuals.Raw; 
   residuals_5.LaggedResiduals(i+1) = reg.residuals(end);
end

residuals_5(1:493,:) =[];
residuals_5(247:end,:) = [];

Predicted_SpotFD_ECM5 = [ECM_5.Coefficients.Estimate(1) + ECM_5.Coefficients.Estimate(3)*data.Fut15_firstdiff_lag1(494:739) ...
    + ECM_5.Coefficients.Estimate(4)*data.Fut15_firstdiff_lag2(494:739) + ECM_5.Coefficients.Estimate(5)*data.Fut15_firstdiff_lag3(494:739) ...
    + ECM_5.Coefficients.Estimate(6)*data.Fut15_firstdiff_lag4(494:739) + ECM_5.Coefficients.Estimate(7)*data.Fut15_firstdiff_lag5(494:739) ...
    + ECM_5.Coefficients.Estimate(8)*data.Fut15_firstdiff_lag6(494:739) + ECM_5.Coefficients.Estimate(9)*data.Fut15_firstdiff_lag7(494:739) ...
    + ECM_5.Coefficients.Estimate(10)*data.Fut15_firstdiff_lag8(494:739) + ECM_5.Coefficients.Estimate(11)*data.Fut15_firstdiff_lag9(494:739) ...
    + ECM_5.Coefficients.Estimate(12)*data.Fut15_firstdiff_lag10(494:739)+ ECM_5.Coefficients.Estimate(13)*data.Spot_firstdiff_lag1(494:739) ...
    + ECM_5.Coefficients.Estimate(14)*data.Spot_firstdiff_lag2(494:739) + ECM_5.Coefficients.Estimate(15)*data.Spot_firstdiff_lag3(494:739) ...
    + ECM_5.Coefficients.Estimate(16)*data.Spot_firstdiff_lag4(494:739) + ECM_5.Coefficients.Estimate(17)*data.Spot_firstdiff_lag5(494:739) ...
    + ECM_5.Coefficients.Estimate(18)*data.Spot_firstdiff_lag6(494:739) + ECM_5.Coefficients.Estimate(19)*data.Spot_firstdiff_lag7(494:739) ...
    + ECM_5.Coefficients.Estimate(20)*data.Spot_firstdiff_lag8(494:739) + ECM_5.Coefficients.Estimate(21)*data.Spot_firstdiff_lag9(494:739) ...
    + ECM_5.Coefficients.Estimate(22)*data.Spot_firstdiff_lag10(494:739)+ ECM_5.Coefficients.Estimate(2)*residuals_5.LaggedResiduals];

Predictions_ECM5 = table();
Predictions_ECM5.Date = data.Date(493:739);
Predictions_ECM5.Spot = data.Spot(493:739);
Predictions_ECM5.PredictedSpotFD = zeros(247,1);
Predictions_ECM5.PredictedSpotFD(2:end) = Predicted_SpotFD_ECM5;
Predictions_ECM5.PredictedSpot = lagmatrix(Predictions_ECM5.Spot, 1) + Predictions_ECM5.PredictedSpotFD;
Predictions_ECM5(1,:) = [];
Predictions_ECM5.Spot = exp(Predictions_ECM5.Spot);
Predictions_ECM5.PredictedSpot = exp(Predictions_ECM5.PredictedSpot);
Predictions_ECM5(238:end,:) = [];

MSE_ECM5 = (sum(Predictions_ECM5.PredictedSpot - Predictions_ECM5.Spot)^2)/(237);
RMSE_ECM5 = sqrt(MSE_ECM5)

ChangePred_ECM5 = Predictions_ECM5.PredictedSpot - lagmatrix(Predictions_ECM5.PredictedSpot, 1);
ChangeSpot_ECM5 = Predictions_ECM5.Spot - lagmatrix(Predictions_ECM5.Spot, 1);

SignPred_ECM5 = table();
for i = 1:height(Predictions_ECM5)
    if ChangePred_ECM5(i) >= 0 & ChangeSpot_ECM5(i) >= 0 
        SignPred_ECM5.True(i) = 1;
    else
        SignPred_ECM5.True(i) = 0;
    end
end

for i = 1:height(Predictions_ECM5)
    if ChangePred_ECM5(i) <= 0 & ChangeSpot_ECM5(i) <= 0 
        SignPred_ECM5.True1(i) = 1;
    else
        SignPred_ECM5.True1(i) = 0;
    end
end

PercentCorrect_ECM5 = (sum(SignPred_ECM5.True) + sum(SignPred_ECM5.True1)) / (height(Predictions_ECM5));



% ECM-GARCH: Spot and Fut16
GARCH5 = garch('GARCHlags', 1, 'ARCHlags', 1)
[estMdl, estParamCov, logL] = estimate(GARCH5, ECM5.Spot_firstdiff)
condVar_5 = infer(estMdl, ECM5.Spot_firstdiff);
ECM5.condVol = sqrt(condVar_5);

ECM_GARCH_5 = fitlm(ECM5, ['Spot_firstdiff ~ LaggedResiduals + Fut15_firstdiff_lag1 + Fut15_firstdiff_lag2' ... 
    '+ Fut15_firstdiff_lag3 + Fut15_firstdiff_lag4 + Fut15_firstdiff_lag5 + Fut15_firstdiff_lag6' ...
    '+ Fut15_firstdiff_lag7 + Fut15_firstdiff_lag8 + Fut15_firstdiff_lag9 + Fut15_firstdiff_lag10' ...
    '+ Spot_firstdiff_lag1 + Spot_firstdiff_lag2 + Spot_firstdiff_lag3 + Spot_firstdiff_lag4' ...
    '+ Spot_firstdiff_lag5 + Spot_firstdiff_lag6+ Spot_firstdiff_lag7 + Spot_firstdiff_lag8' ...
    '+ Spot_firstdiff_lag9 + Spot_firstdiff_lag10 + condVol'])

Pred15 = forecast(estMdl, height(residuals_5), 'Y0', ECM5.Spot_firstdiff);
Pred15_Vol = sqrt(Pred15);


Predicted_SpotFD_GARCH5 = [ECM_GARCH_5.Coefficients.Estimate(1) + ECM_GARCH_5.Coefficients.Estimate(3)*data.Fut15_firstdiff_lag1(494:739) ...
    + ECM_GARCH_5.Coefficients.Estimate(4)*data.Fut15_firstdiff_lag2(494:739) + ECM_GARCH_5.Coefficients.Estimate(5)*data.Fut15_firstdiff_lag3(494:739) ...
    + ECM_GARCH_5.Coefficients.Estimate(6)*data.Fut15_firstdiff_lag4(494:739) + ECM_GARCH_5.Coefficients.Estimate(7)*data.Fut15_firstdiff_lag5(494:739) ...
    + ECM_GARCH_5.Coefficients.Estimate(8)*data.Fut15_firstdiff_lag6(494:739) + ECM_GARCH_5.Coefficients.Estimate(9)*data.Fut15_firstdiff_lag7(494:739) ...
    + ECM_GARCH_5.Coefficients.Estimate(10)*data.Fut15_firstdiff_lag8(494:739) + ECM_GARCH_5.Coefficients.Estimate(11)*data.Fut15_firstdiff_lag9(494:739) ...
    + ECM_GARCH_5.Coefficients.Estimate(12)*data.Fut15_firstdiff_lag10(494:739)+ ECM_GARCH_5.Coefficients.Estimate(13)*data.Spot_firstdiff_lag1(494:739) ...
    + ECM_GARCH_5.Coefficients.Estimate(14)*data.Spot_firstdiff_lag2(494:739) + ECM_GARCH_5.Coefficients.Estimate(15)*data.Spot_firstdiff_lag3(494:739) ...
    + ECM_GARCH_5.Coefficients.Estimate(16)*data.Spot_firstdiff_lag4(494:739) + ECM_GARCH_5.Coefficients.Estimate(17)*data.Spot_firstdiff_lag5(494:739) ...
    + ECM_GARCH_5.Coefficients.Estimate(18)*data.Spot_firstdiff_lag6(494:739) + ECM_GARCH_5.Coefficients.Estimate(19)*data.Spot_firstdiff_lag7(494:739) ...
    + ECM_GARCH_5.Coefficients.Estimate(20)*data.Spot_firstdiff_lag8(494:739) + ECM_GARCH_5.Coefficients.Estimate(21)*data.Spot_firstdiff_lag9(494:739) ...
    + ECM_GARCH_5.Coefficients.Estimate(22)*data.Spot_firstdiff_lag10(494:739) + ECM_GARCH_5.Coefficients.Estimate(23)*Pred15_Vol + ...
    + ECM_GARCH_5.Coefficients.Estimate(2)*residuals_5.LaggedResiduals];

Predictions_GARCH5 = table();
Predictions_GARCH5.Date = data.Date(493:739);
Predictions_GARCH5.Spot = data.Spot(493:739);
Predictions_GARCH5.PredictedSpotFD = zeros(247,1);
Predictions_GARCH5.PredictedSpotFD(2:end) = Predicted_SpotFD_GARCH5;
Predictions_GARCH5.PredictedSpot = lagmatrix(Predictions_GARCH5.Spot, 1) + Predictions_GARCH5.PredictedSpotFD;
Predictions_GARCH5(1,:) = [];
Predictions_GARCH5.Spot = exp(Predictions_GARCH5.Spot);
Predictions_GARCH5.PredictedSpot = exp(Predictions_GARCH5.PredictedSpot);
Predictions_GARCH5(238:end,:) = [];

MSE_GARCH5 = (sum(Predictions_GARCH5.PredictedSpot - Predictions_GARCH5.Spot)^2)/(237);
RMSE_GARCH5 = sqrt(MSE_GARCH5)

ChangePred_GARCH5 = Predictions_GARCH5.PredictedSpot - lagmatrix(Predictions_GARCH5.PredictedSpot, 1);
ChangeSpot_GARCH5 = Predictions_GARCH5.Spot - lagmatrix(Predictions_GARCH5.Spot, 1);

SignPred_GARCH5 = table();
for i = 1:height(Predictions_GARCH5)
    if ChangePred_GARCH5(i) >= 0 & ChangeSpot_GARCH5(i) >= 0 
        SignPred_GARCH5.True(i) = 1;
    else
        SignPred_GARCH5.True(i) = 0;
    end
end

for i = 1:height(Predictions_GARCH5)
    if ChangePred_GARCH5(i) <= 0 & ChangeSpot_GARCH5(i) <= 0 
        SignPred_GARCH5.True1(i) = 1;
    else
        SignPred_GARCH5.True1(i) = 0;
    end
end
PercentCorrect_GARCH5 = (sum(SignPred_GARCH5.True) + sum(SignPred_GARCH5.True1)) / (height(Predictions_GARCH5));


% Finding residuals between Spot and lagged Fut14:
reg6 = table();
reg6.Date = data.Date(1:246);
reg6.Spot = data.Spot(1:246);
reg6.LagFut14 = lagmatrix(data.Fut14(1:246),1);
reg6(1,:) = [];
model6 = fitlm(reg6, 'Spot ~ LagFut14')
reg6.residuals6 = model6.Residuals.Raw;

% Defining ECM variables - Spot and Fut14:
ECM6 = table();
ECM6.Date = reg6.Date;
ECM6.Spot_firstdiff = data.Spot_firstdiff(2:246);
ECM6.LaggedResiduals = lagmatrix(reg6.residuals6,1);
%Lag variables - Fut:
ECM6.Fut14_firstdiff_lag1 = reg6.LagFut14 - lagmatrix(reg6.LagFut14,1);
ECM6.Fut14_firstdiff_lag2 = lagmatrix(ECM6.Fut14_firstdiff_lag1,1);
ECM6.Fut14_firstdiff_lag3 = lagmatrix(ECM6.Fut14_firstdiff_lag1,2);
ECM6.Fut14_firstdiff_lag4 = lagmatrix(ECM6.Fut14_firstdiff_lag1,3);
ECM6.Fut14_firstdiff_lag5 = lagmatrix(ECM6.Fut14_firstdiff_lag1,4);
ECM6.Fut14_firstdiff_lag6 = lagmatrix(ECM6.Fut14_firstdiff_lag1,5);
ECM6.Fut14_firstdiff_lag7 = lagmatrix(ECM6.Fut14_firstdiff_lag1,6);


%Lag variables spot:
ECM6.Spot_firstdiff_lag1 = lagmatrix(ECM6.Spot_firstdiff,1);
ECM6.Spot_firstdiff_lag2 = lagmatrix(ECM6.Spot_firstdiff,2);
ECM6.Spot_firstdiff_lag3 = lagmatrix(ECM6.Spot_firstdiff,3);
ECM6.Spot_firstdiff_lag4 = lagmatrix(ECM6.Spot_firstdiff,4);
ECM6.Spot_firstdiff_lag5 = lagmatrix(ECM6.Spot_firstdiff,5);
ECM6.Spot_firstdiff_lag6 = lagmatrix(ECM6.Spot_firstdiff,6);
ECM6.Spot_firstdiff_lag7 = lagmatrix(ECM6.Spot_firstdiff,7);


ECM_6 = fitlm(ECM6, ['Spot_firstdiff ~ LaggedResiduals + Fut14_firstdiff_lag1 + Fut14_firstdiff_lag2' ... 
    '+ Fut14_firstdiff_lag3 + Fut14_firstdiff_lag4 + Fut14_firstdiff_lag5 + Fut14_firstdiff_lag6' ...
    '+ Fut14_firstdiff_lag7 + Spot_firstdiff_lag1 + Spot_firstdiff_lag2' ...
    '+ Spot_firstdiff_lag3 + Spot_firstdiff_lag4 + Spot_firstdiff_lag5 + Spot_firstdiff_lag6' ...
    '+ Spot_firstdiff_lag7'])

residuals_6 = table();
residuals_6.Date = data.Date;
residuals_6.LaggedResiduals = zeros(1709,1);

for i = 246:493
   reg = table();
   reg.Date = data.Date(1:i);
   reg.Spot = data.Spot(1:i);
   reg.LagFut14 = lagmatrix(data.Fut14(1:i),1);
   reg(1,:) = [];
   model6 = fitlm(reg, 'Spot ~ LagFut14');
   reg.residuals = model6.Residuals.Raw; 
   residuals_6.LaggedResiduals(i+1) = reg.residuals(end);
end

residuals_6(1:246,:) =[];
residuals_6(248:end,:) = [];

Predicted_SpotFD_ECM6 = [ECM_6.Coefficients.Estimate(1) + ECM_6.Coefficients.Estimate(3)*data.Fut14_firstdiff_lag1(247:493) ...
    + ECM_6.Coefficients.Estimate(4)*data.Fut14_firstdiff_lag2(247:493) + ECM_6.Coefficients.Estimate(5)*data.Fut14_firstdiff_lag3(247:493) ...
    + ECM_6.Coefficients.Estimate(6)*data.Fut14_firstdiff_lag4(247:493) + ECM_6.Coefficients.Estimate(7)*data.Fut14_firstdiff_lag5(247:493) ...
    + ECM_6.Coefficients.Estimate(8)*data.Fut14_firstdiff_lag6(247:493) + ECM_6.Coefficients.Estimate(9)*data.Fut14_firstdiff_lag7(247:493) ...
    + ECM_6.Coefficients.Estimate(10)*data.Spot_firstdiff_lag1(247:493) + ECM_6.Coefficients.Estimate(11)*data.Spot_firstdiff_lag2(247:493) ...
    + ECM_6.Coefficients.Estimate(12)*data.Spot_firstdiff_lag3(247:493) + ECM_6.Coefficients.Estimate(13)*data.Spot_firstdiff_lag4(247:493) ...
    + ECM_6.Coefficients.Estimate(14)*data.Spot_firstdiff_lag5(247:493) + ECM_6.Coefficients.Estimate(15)*data.Spot_firstdiff_lag6(247:493) ...
    + ECM_6.Coefficients.Estimate(16)*data.Spot_firstdiff_lag7(247:493) + ECM_6.Coefficients.Estimate(2)*residuals_6.LaggedResiduals];

Predictions_ECM6 = table();
Predictions_ECM6.Date = data.Date(246:493);
Predictions_ECM6.Spot = data.Spot(246:493);
Predictions_ECM6.PredictedSpotFD = zeros(248,1);
Predictions_ECM6.PredictedSpotFD(2:end) = Predicted_SpotFD_ECM6;
Predictions_ECM6.PredictedSpot = lagmatrix(Predictions_ECM6.Spot, 1) + Predictions_ECM6.PredictedSpotFD;
Predictions_ECM6(1,:) = [];
Predictions_ECM6.Spot = exp(Predictions_ECM6.Spot);
Predictions_ECM6.PredictedSpot = exp(Predictions_ECM6.PredictedSpot);
Predictions_ECM6(241:end,:) = [];

MSE_ECM6 = (sum(Predictions_ECM6.PredictedSpot - Predictions_ECM6.Spot)^2)/(240);
RMSE_ECM6 = sqrt(MSE_ECM6)

ChangePred_ECM6 = Predictions_ECM6.PredictedSpot - lagmatrix(Predictions_ECM6.PredictedSpot, 1);
ChangeSpot_ECM6 = Predictions_ECM6.Spot - lagmatrix(Predictions_ECM6.Spot, 1);

SignPred_ECM6 = table();
for i = 1:height(Predictions_ECM6)
    if ChangePred_ECM6(i) >= 0 & ChangeSpot_ECM6(i) >= 0 
        SignPred_ECM6.True(i) = 1;
    else
        SignPred_ECM6.True(i) = 0;
    end
end

for i = 1:height(Predictions_ECM6)
    if ChangePred_ECM6(i) <= 0 & ChangeSpot_ECM6(i) <= 0 
        SignPred_ECM6.True1(i) = 1;
    else
        SignPred_ECM6.True1(i) = 0;
    end
end

PercentCorrect_ECM6 = (sum(SignPred_ECM6.True) + sum(SignPred_ECM6.True1)) / (height(Predictions_ECM6));



% ECM-GARCH: Spot and Fut14:
GARCH6 = garch('GARCHlags', 1, 'ARCHlags', 1)
[estMdl, estParamCov, logL] = estimate(GARCH6, ECM6.Spot_firstdiff)
condVar_6 = infer(estMdl, ECM6.Spot_firstdiff);
ECM6.condVol = sqrt(condVar_6);

ECM_GARCH_6 = fitlm(ECM6, ['Spot_firstdiff ~ LaggedResiduals + Fut14_firstdiff_lag1 + Fut14_firstdiff_lag2' ... 
    '+ Fut14_firstdiff_lag3 + Fut14_firstdiff_lag4 + Fut14_firstdiff_lag5 + Fut14_firstdiff_lag6' ...
    '+ Fut14_firstdiff_lag7 + Spot_firstdiff_lag1 + Spot_firstdiff_lag2 + Spot_firstdiff_lag3' ...
    '+ Spot_firstdiff_lag4 + Spot_firstdiff_lag5 + Spot_firstdiff_lag6+ Spot_firstdiff_lag7' ...
    '+ condVol'])

Pred14 = forecast(estMdl, height(residuals_6), 'Y0', ECM6.Spot_firstdiff);
Pred14_Vol = sqrt(Pred14);


Predicted_SpotFD_GARCH6 = [ECM_GARCH_6.Coefficients.Estimate(1) + ECM_GARCH_6.Coefficients.Estimate(3)*data.Fut14_firstdiff_lag1(247:493) ...
    + ECM_GARCH_6.Coefficients.Estimate(4)*data.Fut14_firstdiff_lag2(247:493) + ECM_GARCH_6.Coefficients.Estimate(5)*data.Fut14_firstdiff_lag3(247:493) ...
    + ECM_GARCH_6.Coefficients.Estimate(6)*data.Fut14_firstdiff_lag4(247:493) + ECM_GARCH_6.Coefficients.Estimate(7)*data.Fut14_firstdiff_lag5(247:493) ...
    + ECM_GARCH_6.Coefficients.Estimate(8)*data.Fut14_firstdiff_lag6(247:493) + ECM_GARCH_6.Coefficients.Estimate(9)*data.Fut14_firstdiff_lag7(247:493) ...
    + ECM_GARCH_6.Coefficients.Estimate(10)*data.Spot_firstdiff_lag1(247:493) + ECM_GARCH_6.Coefficients.Estimate(11)*data.Spot_firstdiff_lag2(247:493) ... 
    + ECM_GARCH_6.Coefficients.Estimate(12)*data.Spot_firstdiff_lag3(247:493) + ECM_GARCH_6.Coefficients.Estimate(13)*data.Spot_firstdiff_lag4(247:493) ...
    + ECM_GARCH_6.Coefficients.Estimate(14)*data.Spot_firstdiff_lag5(247:493) + ECM_GARCH_6.Coefficients.Estimate(15)*data.Spot_firstdiff_lag6(247:493) ...
    + ECM_GARCH_6.Coefficients.Estimate(16)*data.Spot_firstdiff_lag7(247:493) + ECM_GARCH_6.Coefficients.Estimate(17)*Pred14_Vol ...
    + ECM_GARCH_6.Coefficients.Estimate(2)*residuals_6.LaggedResiduals];

Predictions_GARCH6 = table();
Predictions_GARCH6.Date = data.Date(246:493);
Predictions_GARCH6.Spot = data.Spot(246:493);
Predictions_GARCH6.PredictedSpotFD = zeros(248,1);
Predictions_GARCH6.PredictedSpotFD(2:end) = Predicted_SpotFD_GARCH6;
Predictions_GARCH6.PredictedSpot = lagmatrix(Predictions_GARCH6.Spot, 1) + Predictions_GARCH6.PredictedSpotFD;
Predictions_GARCH6(1,:) = [];
Predictions_GARCH6.Spot = exp(Predictions_GARCH6.Spot);
Predictions_GARCH6.PredictedSpot = exp(Predictions_GARCH6.PredictedSpot);
Predictions_GARCH6(241:end,:) = [];

MSE_GARCH6 = (sum(Predictions_GARCH6.PredictedSpot - Predictions_GARCH6.Spot)^2)/(240);
RMSE_GARCH6 = sqrt(MSE_GARCH6)

ChangePred_GARCH6 = Predictions_GARCH6.PredictedSpot - lagmatrix(Predictions_GARCH6.PredictedSpot, 1);
ChangeSpot_GARCH6 = Predictions_GARCH6.Spot - lagmatrix(Predictions_GARCH6.Spot, 1);

SignPred_GARCH6 = table();
for i = 1:height(Predictions_GARCH6)
    if ChangePred_GARCH6(i) >= 0 & ChangeSpot_GARCH6(i) >= 0 
        SignPred_GARCH6.True(i) = 1;
    else
        SignPred_GARCH6.True(i) = 0;
    end
end

for i = 1:height(Predictions_GARCH6)
    if ChangePred_GARCH6(i) <= 0 & ChangeSpot_GARCH6(i) <= 0 
        SignPred_GARCH6.True1(i) = 1;
    else
        SignPred_GARCH6.True1(i) = 0;
    end
end
PercentCorrect_GARCH6 = (sum(SignPred_GARCH6.True) + sum(SignPred_GARCH6.True1)) / (height(Predictions_GARCH6));

% Plots of all models:
figure;
subplot(2,2,1)
plot(Predictions_ECM1.Date, Predictions_ECM1.Spot)
hold on
plot(Predictions_ECM1.Date, Predictions_ECM1.PredictedSpot)
legend('Actual spot', 'Predicted spot', 'Location', 'NorthEast')
title('Predictions by ECM with Dec19 Futures contract')
xlabel('Year')
ylabel('Price in EUR')
box('off')
hold off

subplot(2,2,2)
plot(Predictions_GARCH1.Date, Predictions_GARCH1.Spot)
hold on
plot(Predictions_GARCH1.Date, Predictions_GARCH1.PredictedSpot)
legend('Actual spot', 'Predicted spot', 'Location', 'NorthEast')
title('Predictions by ECM-GARCH with Dec19 Futures contract')
xlabel('Year')
ylabel('Price in EUR')
box('off')
hold off

subplot(2,2,3)
plot(Predictions_ECM2.Date, Predictions_ECM2.Spot)
hold on
plot(Predictions_ECM2.Date, Predictions_ECM2.PredictedSpot)
legend('Actual spot', 'Predicted spot', 'Location', 'NorthEast')
title('Predictions by ECM with Dec18 Futures contract')
xlabel('Year')
ylabel('Price in EUR')
box('off')
hold off

subplot(2,2,4)
plot(Predictions_GARCH2.Date, Predictions_GARCH2.Spot)
hold on
plot(Predictions_GARCH2.Date, Predictions_GARCH2.PredictedSpot)
legend('Actual spot', 'Predicted spot', 'Location', 'NorthEast')
title('Predictions by ECM-GARCH with Dec18 Futures contract')
xlabel('Year')
ylabel('Price in EUR')
box('off')
hold off

figure;
subplot(2,2,1)
plot(Predictions_ECM3.Date, Predictions_ECM3.Spot)
hold on
plot(Predictions_ECM3.Date, Predictions_ECM3.PredictedSpot)
legend('Actual spot', 'Predicted spot', 'Location', 'SouthEast')
title('Predictions by ECM with Dec17 Futures contract')
xlabel('Year')
ylabel('Price in EUR')
box('off')
hold off

subplot(2,2,2)
plot(Predictions_GARCH3.Date, Predictions_GARCH3.Spot)
hold on
plot(Predictions_GARCH3.Date, Predictions_GARCH3.PredictedSpot)
legend('Actual spot', 'Predicted spot', 'Location', 'SouthEast')
title('Predictions by ECM-GARCH with Dec17 Futures contract')
xlabel('Year')
ylabel('Price in EUR')
box('off')
hold off

subplot(2,2,3)
plot(Predictions_ECM4.Date, Predictions_ECM4.Spot)
hold on
plot(Predictions_ECM4.Date, Predictions_ECM4.PredictedSpot)
legend('Actual spot', 'Predicted spot', 'Location', 'NorthEast')
title('Predictions by ECM with Dec16 Futures contract')
xlabel('Year')
ylabel('Price in EUR')
box('off')
hold off

subplot(2,2,4)
plot(Predictions_GARCH4.Date, Predictions_GARCH4.Spot)
hold on
plot(Predictions_GARCH4.Date, Predictions_GARCH4.PredictedSpot)
legend('Actual spot', 'Predicted spot', 'Location', 'NorthEast')
title('Predictions by ECM-GARCH with Dec16 Futures contract')
xlabel('Year')
ylabel('Price in EUR')
box('off')
hold off

figure;
subplot(2,2,1)
plot(Predictions_ECM5.Date, Predictions_ECM5.Spot)
hold on
plot(Predictions_ECM5.Date, Predictions_ECM5.PredictedSpot)
legend('Actual spot', 'Predicted spot', 'Location', 'SouthEast')
title('Predictions by ECM with Dec15 Futures contract')
xlabel('Year')
ylabel('Price in EUR')
box('off')
hold off

subplot(2,2,2)
plot(Predictions_GARCH5.Date, Predictions_GARCH5.Spot)
hold on
plot(Predictions_GARCH5.Date, Predictions_GARCH5.PredictedSpot)
legend('Actual spot', 'Predicted spot', 'Location', 'SouthEast')
title('Predictions by ECM-GARCH with Dec15 Futures contract')
xlabel('Year')
ylabel('Price in EUR')
box('off')
hold off

subplot(2,2,3)
plot(Predictions_ECM6.Date, Predictions_ECM6.Spot)
hold on
plot(Predictions_ECM6.Date, Predictions_ECM6.PredictedSpot)
legend('Actual spot', 'Predicted spot', 'Location', 'SouthEast')
title('Predictions by ECM with Dec14 Futures contract')
xlabel('Year')
ylabel('Price in EUR')
box('off')
hold off

subplot(2,2,4)
plot(Predictions_GARCH6.Date, Predictions_GARCH6.Spot)
hold on
plot(Predictions_GARCH6.Date, Predictions_GARCH6.PredictedSpot)
legend('Actual spot', 'Predicted spot', 'Location', 'SouthEast')
title('Predictions by ECM-GARCH with Dec14 Futures contract')
xlabel('Year')
ylabel('Price in EUR')
box('off')
hold off

%% Part 6 - Link to fundamental variables:
clear;
clc;

% Import data and create variables:
data = readtable('Spot+DecFutures.xlsx');
data=data(~any(ismissing(data.Spot),2),:);
data(584:585,:) = [];
data = table2timetable(data);
data.Spot = log(data.Spot);
data.Oil = log(data.Oil);
data.Electricity = log(data.Electricity);
data.Elec = data.Electricity;
data.Coal = log(data.Coal);
data.DAX = log(data.DAX);
data.NaturalGAS = log(data.NaturalGAS);
data.Ngas = data.NaturalGAS;
data(1:556,:) = [];


% First differences:
data.Spot_fd = data.Spot - lagmatrix(data.Spot,1);
data.Oil_fd = data.Oil - lagmatrix(data.Oil,1);
data.Elec_fd = data.Elec - lagmatrix(data.Elec,1);
data.Coal_fd = data.Coal - lagmatrix(data.Coal,1);
data.DAX_fd = data.DAX - lagmatrix(data.DAX,1);
data.Ngas_fd = data.Ngas - lagmatrix(data.Ngas,1);


% Determining lag length for VAR (full model):
reg_VAR1 = table();
reg_VAR1.Spot_fd = data.Spot_fd;
reg_VAR1.Oil_fd = data.Oil_fd;
reg_VAR1.Elec_fd = data.Elec_fd;
reg_VAR1.Coal_fd = data.Coal_fd;
reg_VAR1.DAX_fd = data.DAX_fd;
reg_VAR1.Ngas_fd = data.Ngas_fd;
reg_VAR1(1707:end,:) = [];
reg_VAR1=reg_VAR1(~any(ismissing(reg_VAR1.Elec_fd),2),:);
reg_VAR1=reg_VAR1(~any(ismissing(reg_VAR1.DAX_fd),2),:);
reg_VAR1=reg_VAR1(~any(ismissing(reg_VAR1.Ngas_fd),2),:);
reg_VAR1=reg_VAR1(~any(ismissing(reg_VAR1.Coal_fd),2),:);
reg_VAR1 = reg_VAR1{:,:};
[laglength, AIC, logL] = VARlag(reg_VAR1,12)
%Optimal = 2


% Simple regressions:
data1 = table();
data1.Spot_fd = data.Spot_fd;
data1.Oil_fd = data.Oil_fd;
data1.Elec_fd = data.Elec_fd;
data1.Coal_fd = data.Coal_fd;
data1.DAX_fd = data.DAX_fd;
data1.Ngas_fd = data.Ngas_fd;
OLS_1 = fitlm(data1, 'Spot_fd ~ Oil_fd + Elec_fd + Coal_fd + DAX_fd + Ngas_fd')


% Spot carbon and oil:
data2 = table();
data2.Spot_fd = data.Spot_fd;
data2.Oil_fd = data.Oil_fd;
OLS_2 = fitlm(data2, 'Spot_fd ~ Oil_fd')


% Spot carbon and elec:
data3 = table();
data3.Spot_fd = data.Spot_fd;
data3.Elec_fd = data.Elec_fd;
OLS_3 = fitlm(data3, 'Spot_fd ~ Elec_fd')

% Spot carbon and Coal:
data4 = table();
data4.Spot_fd = data.Spot_fd;
data4.Coal_fd = data.Coal_fd;
OLS_4 = fitlm(data4, 'Spot_fd ~ Coal_fd')

% Spot carbon and DAX:
data5 = table();
data5.Spot_fd = data.Spot_fd;
data5.DAX_fd = data.DAX_fd;
OLS_5 = fitlm(data5, 'Spot_fd ~ DAX_fd')

% Spot and Ngas:
data6 = table();
data6.Spot_fd = data.Spot_fd;
data6.Ngas_fd = data.Ngas_fd;
OLS_6 = fitlm(data6, 'Spot_fd ~ Ngas_fd')

% Second analysis - Full VAR with impulse responses:
table1 = data(:, {'Spot_fd', 'Oil_fd', 'Elec_fd', 'Coal_fd', 'DAX_fd', 'Ngas_fd'});
k = 2;
g = size(table1,2);
model = estimate( varm( g , k ), table1.Variables );
results = summarize(model);

% Extracting summary statistics to excel:
writetable(results.Table, 'VARtest2.xlsx', 'Sheet', 1, 'Range', 'B1');

T = 20;
Y = armairf(model.AR, [], 'NumObs', T);

% plot
figure('WindowStyle', 'docked');
subplot(2,3,1); plot((1:T)', Y(:,1,1), 'LineWidth', 1); title('Response of Spot to Spot'); grid('on');
subplot(2,3,2); plot((1:T)', Y(:,1,2), 'LineWidth', 1); ylim([-0.06, 0.06]); title('Response of Spot to Oil'); grid('on');
subplot(2,3,3); plot((1:T)', Y(:,1,3), 'LineWidth', 1); ylim([-0.01, 0.01]); title('Response of Spot to Elec'); grid('on');
subplot(2,3,4); plot((1:T)', Y(:,1,4), 'LineWidth', 1); title('Response of Spot to Coal'); grid('on');
subplot(2,3,5); plot((1:T)', Y(:,1,5), 'LineWidth', 1); ylim([-0.03, 0.07]); title('Response of Spot to DAX'); grid('on');
subplot(2,3,6); plot((1:T)', Y(:,1,6), 'LineWidth', 1); ylim([-0.02, 0.01]); title('Response of Spot to Ngas'); grid('on');
%%