%% Master script


%% Load data
NOKUSDTable = readtable('NOKUSD.xlsx');
 
OilFuturesPricesTable = readtable('OilFuturesPrices.xlsx');
 
%% convert tables into arrays without the date-coloumn
NOKUSD = table2array(NOKUSDTable(:,2));
 
OilFuturesPrices = table2array(OilFuturesPricesTable(:,2:13));

%% Create the termstructur of oil.

TermStructureOil = (OilFuturesPrices(:, 2:end) - OilFuturesPrices(:, 1))./OilFuturesPrices(:, 1);

%% Currency return
CurrencyReturn = (NOKUSD(2:end) - NOKUSD(1:end-1))./NOKUSD(1:end-1);
CurRet = CurrencyReturn*100;

%% DieboldLi factor loadings and betas
TSO = [OilFuturesPrices(:, 1), OilFuturesPrices(:, 1), OilFuturesPrices(:, 1), OilFuturesPrices(:, 1), OilFuturesPrices(:, 1), OilFuturesPrices(:, 1), OilFuturesPrices(:, 1), OilFuturesPrices(:, 1), OilFuturesPrices(:, 1), OilFuturesPrices(:, 1), OilFuturesPrices(:, 1)];
TermStructureOil = (OilFuturesPrices(:, 2:end) - TSO)./TSO;


%% DieboldLi factor loadings and betas
% Need to use matlab version from 2016.
EstimationData = TermStructureOil(2:end,:);

% Explicitly set the time factor lambda
lambda_t = .0609;

% Construct a matrix of the factor loadings
% Tenors associated with data
TimeToMat = [1 2 3 4 5 6 7 8 9 10 11]';
X = [ones(size(TimeToMat)) (1 - exp(-lambda_t*TimeToMat))./(lambda_t*TimeToMat) ...
    ((1 - exp(-lambda_t*TimeToMat))./(lambda_t*TimeToMat) - exp(-lambda_t*TimeToMat))];

% Plot the factor loadings
figure
set(gcf,'color','w');
plot(TimeToMat,X)
title('Factor Loadings for Diebold Li Model with time factor of .0609')
xlabel('Maturity (months)')
ylim([0 1.1])
legend({'Level','Slope','Curvature'},'location','east')




% Preallocate the Betas
Beta = zeros(size(EstimationData,1),3);

% Loop through and fit each end of month yield curve
for jdx = 1:size(EstimationData,1)
    tmpCurveModel = DieboldLi.fitBetasFromYields(jdx,lambda_t*12,daysadd(jdx,30*TimeToMat),EstimationData(jdx,:)');
    Beta(jdx,:) = [tmpCurveModel.Beta1 tmpCurveModel.Beta2 tmpCurveModel.Beta3];
end


%% Plots: USDNOK against First near by of oil and others
figure
set(gcf,'color','w');
plot(NOKUSD,'k')
title('PLOT OF NOK/USD AND FIRST-NEAR-BY CRUDE BRENT OIL FUTURES PRICE');
xlabel('Year');
ylabel('NOK/USD');
hold on
yyaxis right
plot(OilFuturesPrices(:,1),'--')
ylabel('FIRST NEAR BY OIL FUTURES PRICE IN USD');
legend('NOK/USD','FNB OIL FUTURES PRICE','Location','northwest')
xticks(1:12:242)
xticklabels({2001:2021})
xtickangle(45)
xlim([4 243])


%% Plot the termstructure of oil. 

 
dates = 1:243;
maturities = (1:11)';
values=TermStructureOil;
% Plotting the mesh plot.

figure
set(gcf,'color','w');
mesh(dates,maturities,values');
xlabel('YEAR');
ylabel('BASIS');
zlabel('%-WISE VARIATION');
set (gca, 'XLim', [dates(1), dates(end)])
set (gca, 'YLim', [1 11])
xticks(1:12:242)
xticklabels({2001:2021})
xtickangle(45)
xlim([4 243])


%% Plots: Betas 

figure
set(gcf,'color','w');
subplot(3,1,1);
plot(Betas(1:end,1), 'LineWidth', 1)
title('Estimated Level Factor');
xticks(1:12:243)
xticklabels({2001:2021})
xtickangle(45)
xlim([4 243])
set(gca,'box','off')

subplot(3,1,2);
plot(Betas(1:end,2), 'LineWidth', 1)
title('Estimated Slope Factor');
xticks(1:12:243)
xticklabels({2001:2021})
xtickangle(45)
xlim([4 243])
set(gca,'box','off')

subplot(3,1,3);
plot(Betas(1:end,3), 'LineWidth', 1)
title('Estimated Curvature Factor');
xticks(1:12:243)
xticklabels({2001:2021})
xtickangle(45)
xlim([4 243])
set(gca,'box','off')


%% Benchmark Random walk without drift

y = CurRet;

yhat_RW = zeros(length(y),1);
yact_RW = y;
err_RW = yact_RW(61:end) - yhat_RW(61:end);
SFE_RW = err_RW.^2;
AFE_RW = abs(err_RW);

MSFE_RW = mean(SFE_RW);
MAFE_RW = mean(AFE_RW);

RMSFE_RW = sqrt(MSFE_RW)

dates = (1:size(yhat_RW-60,1))';
figure;
subplot(3,1,1); plot(dates(1:end-60),[yact_RW(61:end) yhat_RW(61:end)]); legend('actual','forecast'); title('Constant Forecast')
subplot(3,1,2); plot(dates(1:end-60),SFE_RW); title('SFE')
subplot(3,1,3); plot(dates(1:end-60),AFE_RW); title('AFE')

%% Random walk with drift rolling
win = 60;
yhat_RWdrift = zeros(length(y)-win,1);
for i = 1:length(y)-win 
    ywin = y(i:i+win-1);
    yhat_RWdrift(i) = sum(ywin)/win;
end


yact_RWdrift = y(61:end);
err_RWdrift = yact_RWdrift - yhat_RWdrift;
SFE_RWdrift = err_RWdrift.^2;
AFE_RWdrift = abs(err_RWdrift);

MSFE_RWdrift = mean(SFE_RWdrift);
MAFE_RWdrift = mean(AFE_RWdrift);

RMSFE_RWdrift = sqrt(MSFE_RWdrift)

dates = (1:size(yhat_RWdrift,1))';
figure;
subplot(3,1,1); plot(dates,[yact_RWdrift yhat_RWdrift]); legend('actual','forecast'); title('Forecast RW with drift')
subplot(3,1,2); plot(dates,SFE_RWdrift); title('SFE RW with drift')
subplot(3,1,3); plot(dates,AFE_RWdrift); title('AFE RW with drift')


%% Random walk with drift expanding
win = 60;
yhat_RWdrift = zeros(length(y)-win,1);
for i = 1:length(y)-win 
    ywin = y(1:i+win-1);
    yhat_RWdrift(i) = sum(ywin)/(win+i-1);
end


yact_RWdrift = y(61:end);
err_RWdrift = yact_RWdrift - yhat_RWdrift;
SFE_RWdrift = err_RWdrift.^2;
AFE_RWdrift = abs(err_RWdrift);

MSFE_RWdrift = mean(SFE_RWdrift);
MAFE_RWdrift = mean(AFE_RWdrift);

RMSFE_RWdrift = sqrt(MSFE_RWdrift)

dates = (1:size(yhat_RWdrift,1))';
figure;
subplot(3,1,1); plot(dates,[yact_RWdrift yhat_RWdrift]); legend('actual','forecast'); title('Forecast RW with drift')
subplot(3,1,2); plot(dates,SFE_RWdrift); title('SFE RW with drift')
subplot(3,1,3); plot(dates,AFE_RWdrift); title('AFE RW with drift')

%% ARX rolling window, remove xt for linear regression forecast
p = 0;
win = 60;
[T,nvar] = size(y);
syact_AR_win = zeros(T-win,1); % store actual values by row
sfcasts_AR_win = zeros(T-win,1); % store forecasts by row
serr_AR_win = zeros(T-win,1); % store errors by row

for t = 1:T-win
    % get y and x values
    ywin = y(t:t+win-1); % get y for window
    [nobs,nvar] = size(ywin); 
    xt = zeros(nobs-p,p);
    for ii=1:p
        xt(:,(nvar*(ii-1)+1):nvar*ii)=ywin(ii:nobs-p+ii-1,:);
    end     
    wt = Betas(2+t-1:win+t-p,1:end)
    x=[ones(nobs-p,1) xt wt]; 
    % remove inital lags from y
    yshort_AR_win = ywin(p+1:end,:);
    % compute beta
    betahat_AR_win = (x'*x)\(x'*yshort_AR_win);
    % compute forecast
    yhat_AR_win = x(end,:)*betahat_AR_win;
    % compute forecast error
    yact = y(t+win);
    err_AR_win = yact - yhat_AR_win; 
    
    % Store
    sfcasts_AR_win(t,:) = yhat_AR_win;
    syact(t,:) = yact;
    serr_AR_win(t,:) = err_AR_win; 
end
SFE_AR_win = serr_AR_win.^2;
AFE_AR_win = abs(serr_AR_win);
CSSED_AR_win = cumsum(SFE_RW - SFE_AR_win);

MSFE_AR_win = mean(SFE_AR_win);
MAFE_AR_win = mean(AFE_AR_win);

RMSFE_AR_win = sqrt(MSFE_AR_win)

dates = (1:T-win)';
figure;
subplot(3,1,1); plot(dates,[syact sfcasts_AR_win]); legend('actual','forecast'); title('Rolling Window AR(p) Forecast')
subplot(3,1,2); plot(dates,SFE_AR_win); title('SFE')
subplot(3,1,3); plot(dates,AFE_AR_win); title('AFE')

figure;
plot(CSSED_AR_win)

%% ARX expanding window, remove xt for linear regression forecast
p = 0;
win = 60;
[T,nvar] = size(y);
syact_AR_win = zeros(T-win,1); % store actual values by row
sfcasts_AR_win = zeros(T-win,1); % store forecasts by row
serr_AR_win = zeros(T-win,1); % store errors by row

for t = 1:T-win
    % get y and x values
    ywin = y(1:t+win-1); % get y for window
    [nobs,nvar] = size(ywin); 
    xt = zeros(nobs-p,p);
    for ii=1:p
        xt(:,(nvar*(ii-1)+1):nvar*ii)=ywin(ii:nobs-p+ii-1,:);
    end     
    wt = Betas(2+1-1:win+t-p,1:end)
    x=[ones(nobs-p,1) xt wt]; 
    % remove inital lags from y
    yshort_AR_win = ywin(p+1:end,:);
    % compute beta
    betahat_AR_win = (x'*x)\(x'*yshort_AR_win);
    % compute forecast
    yhat_AR_win = x(end,:)*betahat_AR_win;
    % compute forecast error
    yact = y(t+win);
    err_AR_win = yact - yhat_AR_win; 
    
    % Store
    sfcasts_AR_win(t,:) = yhat_AR_win;
    syact(t,:) = yact;
    serr_AR_win(t,:) = err_AR_win; 
end
SFE_AR_win = serr_AR_win.^2;
AFE_AR_win = abs(serr_AR_win);
CSSED_AR_win = cumsum(SFE_RW - SFE_AR_win);

MSFE_AR_win = mean(SFE_AR_win);
MAFE_AR_win = mean(AFE_AR_win);

RMSFE_AR_win = sqrt(MSFE_AR_win)

dates = (1:T-win)';
figure;
subplot(3,1,1); plot(dates,[syact sfcasts_AR_win]); legend('actual','forecast'); title('Rolling Window AR(p) Forecast')
subplot(3,1,2); plot(dates,SFE_AR_win); title('SFE')
subplot(3,1,3); plot(dates,AFE_AR_win); title('AFE')

figure;
plot(CSSED_AR_win)


%% Rolling Window AR(p)
p = 1;
win = 60;
[T,nvar] = size(y);
syact_AR_win = zeros(T-win,1); % store actual values by row
sfcasts_AR_win = zeros(T-win,1); % store forecasts by row
serr_AR_win = zeros(T-win,1); % store errors by row

for t = 1:T-win
    % get y and x values
    ywin = y(t:t+win-1); % get y for window
    [nobs,nvar] = size(ywin); 
    xt = zeros(nobs-p,p);
    for ii=1:p
        xt(:,(nvar*(ii-1)+1):nvar*ii)=ywin(ii:nobs-p+ii-1,:);
    end     
    x=[ones(nobs-p,1) xt]; 
    % remove inital lags from y
    yshort_AR_win = ywin(p+1:end,:);
    % compute beta
    betahat_AR_win = (x'*x)\(x'*yshort_AR_win);
    % compute forecast
    yhat_AR_win = x(end,:)*betahat_AR_win;
    % compute forecast error
    yact = y(t+win);
    err_AR_win = yact - yhat_AR_win; 
    
    % Store
    sfcasts_AR_win(t,:) = yhat_AR_win;
    syact(t,:) = yact;
    serr_AR_win(t,:) = err_AR_win; 
end
SFE_AR_win = serr_AR_win.^2;
AFE_AR_win = abs(serr_AR_win);
CSSED_AR_win = cumsum(SFE_RW - SFE_AR_win);

MSFE_AR_win = mean(SFE_AR_win);
MAFE_AR_win = mean(AFE_AR_win);

RMSFE_AR_win = sqrt(MSFE_AR_win)

dates = (1:T-win)';
figure;
subplot(3,1,1); plot(dates,[syact sfcasts_AR_win]); legend('actual','forecast'); title('Rolling Window AR(p) Forecast')
subplot(3,1,2); plot(dates,SFE_AR_win); title('SFE')
subplot(3,1,3); plot(dates,AFE_AR_win); title('AFE')

figure;
plot(CSSED_AR_win)


%% Expanding Window AR(p)
p = 10;
win = 60;
[T,nvar] = size(y);
syact_AR_win = zeros(T-win,1); % store actual values by row
sfcasts_AR_win = zeros(T-win,1); % store forecasts by row
serr_AR_win = zeros(T-win,1); % store errors by row

for t = 1:T-win
    % get y and x values
    ywin = y(1:t+win-1); % get y for window
    [nobs,nvar] = size(ywin); 
    xt = zeros(nobs-p,p);
    for ii=1:p
        xt(:,(nvar*(ii-1)+1):nvar*ii)=ywin(ii:nobs-p+ii-1,:);
    end     
    x=[ones(nobs-p,1) xt]; 
    % remove inital lags from y
    yshort_AR_win = ywin(p+1:end,:);
    % compute beta
    betahat_AR_win = (x'*x)\(x'*yshort_AR_win);
    % compute forecast
    yhat_AR_win = x(end,:)*betahat_AR_win;
    % compute forecast error
    yact = y(t+win);
    err_AR_win = yact - yhat_AR_win; 
    
    % Store
    sfcasts_AR_win(t,:) = yhat_AR_win;
    syact(t,:) = yact;
    serr_AR_win(t,:) = err_AR_win; 
end
SFE_AR_win = serr_AR_win.^2;
AFE_AR_win = abs(serr_AR_win);
CSSED_AR_win = cumsum(SFE_RW - SFE_AR_win);

MSFE_AR_win = mean(SFE_AR_win);
MAFE_AR_win = mean(AFE_AR_win);

RMSFE_AR_win = sqrt(MSFE_AR_win)

dates = (1:T-win)';
figure;
subplot(3,1,1); plot(dates,[syact sfcasts_AR_win]); legend('actual','forecast'); title('Rolling Window AR(p) Forecast')
subplot(3,1,2); plot(dates,SFE_AR_win); title('SFE')
subplot(3,1,3); plot(dates,AFE_AR_win); title('AFE')

figure;
plot(CSSED_AR_win)


%% LSTM currencyreturn
data = CurRet';



YPred = [];


for i = 1:182

numTimeStepsTrain = 59+i

dataTrain = data(1:numTimeStepsTrain+1);
dataTest = data(numTimeStepsTrain+1);


%
mu = mean(dataTrain);
sig = std(dataTrain);

dataTrainStandardized = (dataTrain - mu) / sig;


%
XTrain = dataTrainStandardized(1:end-1);
YTrain = dataTrainStandardized(2:end);


%
numFeatures = 1;
numResponses = 1;
numHiddenUnits = 200;

layers = [ ...
    sequenceInputLayer(numFeatures)
    lstmLayer(numHiddenUnits)
    fullyConnectedLayer(numResponses)
    regressionLayer];


%
options = trainingOptions('adam', ...
    'MaxEpochs',60, ...
    'GradientThreshold',1, ...
    'InitialLearnRate',0.005, ...
    'LearnRateSchedule','piecewise', ...
    'LearnRateDropPeriod',125, ...
    'LearnRateDropFactor',0.2, ...
    'Verbose',0, ...
    'Plots','training-progress');



net = trainNetwork(XTrain,YTrain,layers,options);


%
dataTestStandardized = (dataTest - mu) / sig;
XTest = dataTestStandardized;

% Reset network and give X values
net = resetState(net);
net = predictAndUpdateState(net,XTrain);

close all force

%
[net,YPred(i)] = predictAndUpdateState(net,XTest,'ExecutionEnvironment','cpu');
YPred(i) = sig*YPred(i) + mu;
end


%%
YTest = data(61:end);
err_LSTM = YTest-YPred;
SFE_LSTM_Currency = (YTest-YPred).^2;
rmse = sqrt(mean((YTest-YPred).^2))
CSSED_LSTM_Currency = cumsum(SFE_RW - SFE_LSTM_Currency');

%%
figure
subplot(2,1,1)
plot(YTest)
hold on
plot(YPred,'.-')
hold off
legend(["Observed" "Predicted"])
ylabel("Return")
title("Forecast with Updates")

subplot(2,1,2)
stem(YPred - YTest)
xlabel("Month")
ylabel("Error")
title("RMSE = " + rmse)

figure;
plot(CSSED_LSTM_Currency)



%% Expanding LSTM Betas

% Load data
dataY = CurRet';
dataX = Betas(2:end,:)';


% Separate training- and test-data.


YPred = [];


for i = 1:182

numTimeStepsTrain = 59+i

%
dataTrainX = dataX(1:numTimeStepsTrain+1);
dataTestX = dataX(numTimeStepsTrain+1);

dataTrainY = dataY(1:numTimeStepsTrain+1);
dataTestY = dataY(numTimeStepsTrain+1);

%
muX = mean(dataTrainX);
sigX = std(dataTrainX);

muY = mean(dataTrainY);
sigY = std(dataTrainY);

%
dataTrainStandardizedX = (dataTrainX - muX) / sigX;

dataTrainStandardizedY = (dataTrainY - muY) / sigY;

%
XTrain = dataTrainStandardizedX(1:end-1);
YTrain = dataTrainStandardizedY(2:end);


%
numFeatures = 1;
numResponses = 1;
numHiddenUnits = 200;

layers = [ ...
    sequenceInputLayer(numFeatures)
    lstmLayer(numHiddenUnits)
    fullyConnectedLayer(numResponses)
    regressionLayer];


%
options = trainingOptions('adam', ...
    'MaxEpochs',50, ...
    'GradientThreshold',1, ...
    'InitialLearnRate',0.005, ...
    'LearnRateSchedule','piecewise', ...
    'LearnRateDropPeriod',125, ...
    'LearnRateDropFactor',0.2, ...
    'Verbose',0, ...
    'Plots','training-progress');



net = trainNetwork(XTrain,YTrain,layers,options);


%
dataTestStandardizedX = (dataTestX - muX) / sigX;
XTest = dataTestStandardizedX;

% Expanding window forecasting
net = resetState(net);
net = predictAndUpdateState(net,XTrain);

close all force

%
[net,YPred(i)] = predictAndUpdateState(net,XTest,'ExecutionEnvironment','cpu');
YPred(i) = sigY*YPred(i) + muY;
end




%%
YTest = dataY(61:end);
err_LSTM_Betas = YTest-YPred;
SFE_LSTM_Currency = (YTest-YPred).^2;
rmse = sqrt(mean((YTest-YPred).^2))
CSSED_LSTM_Currency = cumsum(SFE_RW - SFE_LSTM_Currency');


%%
figure
subplot(2,1,1)
plot(YTest)
hold on
plot(YPred,'.-')
hold off
legend(["Observed" "Predicted"])
ylabel("Return")
title("Forecast with Updates")

subplot(2,1,2)
stem(YPred - YTest)
xlabel("Month")
ylabel("Error")
title("RMSE = " + rmse)

figure;
plot(CSSED_LSTM_Currency)



%% LSTM BIG DATASET

BigDatasetTable = readtable('BigDataset.xlsx');
 

 
%% convert tables into arrays without the date-coloumn
BigDataset = table2array(BigDatasetTable(:,2:end));

BigDataset = [BigDataset, Betas];


%% Rolling LSTM BigDataset

% Load data
dataY = CurRet';
dataX = BigDataset(2:end,:)';


% Separate training- and test-data.


YPred = [];


for i = 1:182

numTimeStepsTrain = 59+i

%
dataTrainX = dataX(1:numTimeStepsTrain+1);
dataTestX = dataX(numTimeStepsTrain+1);

dataTrainY = dataY(1:numTimeStepsTrain+1);
dataTestY = dataY(numTimeStepsTrain+1);

%
muX = mean(dataTrainX);
sigX = std(dataTrainX);

muY = mean(dataTrainY);
sigY = std(dataTrainY);

%
dataTrainStandardizedX = (dataTrainX - muX) / sigX;

dataTrainStandardizedY = (dataTrainY - muY) / sigY;

%
XTrain = dataTrainStandardizedX(1:end-1);
YTrain = dataTrainStandardizedY(2:end);


%
numFeatures = 1;
numResponses = 1;
numHiddenUnits = 200;

layers = [ ...
    sequenceInputLayer(numFeatures)
    lstmLayer(numHiddenUnits)
    fullyConnectedLayer(numResponses)
    regressionLayer];


%
options = trainingOptions('adam', ...
    'MaxEpochs',50, ...
    'GradientThreshold',1, ...
    'InitialLearnRate',0.005, ...
    'LearnRateSchedule','piecewise', ...
    'LearnRateDropPeriod',125, ...
    'LearnRateDropFactor',0.2, ...
    'Verbose',0, ...
    'Plots','training-progress');



net = trainNetwork(XTrain,YTrain,layers,options);


%
dataTestStandardizedX = (dataTestX - muX) / sigX;
XTest = dataTestStandardizedX;

% Expanding window forecasting
net = resetState(net);
net = predictAndUpdateState(net,XTrain);

close all force

%
[net,YPred(i)] = predictAndUpdateState(net,XTest,'ExecutionEnvironment','cpu');
YPred(i) = sigY*YPred(i) + muY;
end




%%
YTest = dataY(61:end);
err_LSTM_Betas = YTest-YPred;
SFE_LSTM_Currency = (YTest-YPred).^2;
rmse = sqrt(mean((YTest-YPred).^2))
CSSED_LSTM_Currency = cumsum(SFE_RW - SFE_LSTM_Currency');


%%
figure
subplot(2,1,1)
plot(YTest)
hold on
plot(YPred,'.-')
hold off
legend(["Observed" "Predicted"])
ylabel("Return")
title("Forecast with Updates")

subplot(2,1,2)
stem(YPred - YTest)
xlabel("Month")
ylabel("Error")
title("RMSE = " + rmse)

figure;
plot(CSSED_LSTM_Currency)



%% Model combination traditional rolling
yhat_MClrolling = 0.25*yhat_AR1rolling + 0.25*yhat_ARX1rolling + 0.25*yhat_Linearrolling + 0.25*yhat_RWdriftrolling;

err_MCrolling = CurRet(61:end) - yhat_MCrolling;

SFE_AR_win = err_MCrolling.^2;
CSSED_AR_win = cumsum(SFE_RW - SFE_AR_win);

MSFE_AR_win = mean(SFE_AR_win);

RMSFE_AR_win = sqrt(MSFE_AR_win)


figure;
plot(CSSED_AR_win)


%% Model combination traditional expanding
yhat_MCtraditionalexpanding = 0.25*yhat_AR1expanding + 0.25*yhat_ARX1expanding + 0.25*yhat_Linearexpanding + 0.25*yhat_RWdriftexpanding;

err_MCexpanding = CurRet(61:end) - yhat_MCexpanding;

SFE_AR_win = err_MCexpanding.^2;
CSSED_AR_win = cumsum(SFE_RW - SFE_AR_win);

MSFE_AR_win = mean(SFE_AR_win);

RMSFE_AR_win = sqrt(MSFE_AR_win)

figure;
plot(CSSED_AR_win)




%% Model combination betas rolling 
yhat_MCrolling = 0.2*yhat_AR1rolling + 0.2*yhat_ARX1rolling + 0.2*yhat_Linearrolling + 0.2*yhat_LSTMbetasrolling' + 0.2*yhat_RWdriftrolling;

err_MCrolling = CurRet(61:end) - yhat_MCrolling;

SFE_AR_win = err_MCrolling.^2;
CSSED_AR_win = cumsum(SFE_RW - SFE_AR_win);

MSFE_AR_win = mean(SFE_AR_win);

RMSFE_AR_win = sqrt(MSFE_AR_win)


figure;
plot(CSSED_AR_win)


%% Model combination betas expanding
yhat_MCexpanding = 0.2*yhat_AR1expanding + 0.2*yhat_ARX1expanding + 0.2*yhat_Linearexpanding + 0.2*yhat_LSTMbetasexpanding' + 0.2*yhat_RWdriftexpanding;

err_MCexpanding = CurRet(61:end) - yhat_MCexpanding;

SFE_AR_win = err_MCexpanding.^2;
CSSED_AR_win = cumsum(SFE_RW - SFE_AR_win);

MSFE_AR_win = mean(SFE_AR_win);

RMSFE_AR_win = sqrt(MSFE_AR_win)

figure;
plot(CSSED_AR_win)


%% Model combination big dataset rolling
yhat_MCrolling = 0.2*yhat_AR1rolling + 0.2*yhat_ARX1rolling + 0.2*yhat_Linearrolling + 0.2*yhat_LSTMbigdatarolling' + 0.2*yhat_RWdriftrolling;

err_MCrolling = CurRet(61:end) - yhat_MCrolling;

SFE_AR_win = err_MCrolling.^2;
CSSED_AR_win = cumsum(SFE_RW - SFE_AR_win);

MSFE_AR_win = mean(SFE_AR_win);

RMSFE_AR_win = sqrt(MSFE_AR_win)

figure;
plot(CSSED_AR_win)


%% Model combination big dataset expanding
yhat_MCbigdataexpanding = 0.2*yhat_AR1expanding + 0.2*yhat_ARX1expanding + 0.2*yhat_Linearexpanding + 0.2*yhat_LSTMbigdataexpanding' + 0.2*yhat_RWdriftexpanding;

err_MCexpanding = CurRet(61:end) - yhat_MCexpanding;

SFE_AR_win = err_MCexpanding.^2;
CSSED_AR_win = cumsum(SFE_RW - SFE_AR_win);

MSFE_AR_win = mean(SFE_AR_win);

RMSFE_AR_win = sqrt(MSFE_AR_win)

figure;
plot(CSSED_AR_win)



%% Diebold Mariano west test

% Initialization
err1 = err_RWnodrift;
err2 = err_RWdriftrolling;
err3 = err_RWdriftexpanding;
err4 = err_AR1rolling;
err5 = err_AR1expanding;
err6 = err_Linearrolling;
err7 = err_Linearexpanding;
err8 = err_ARX1rolling;
err9 = err_ARX1expanding;
err10 = err_LSTMbetasrolling';
err11 = err_LSTMbetasexpanding';
err12 = err_MCrolling;
err13 = err_MCexpanding;
err14 = err_LSTMcurretrolling';
err15 = err_LSTMcurretexpanding';
err16 = err_LSTMbigdatarolling';
err17 = err_LSTMbigdataexpanding';


e1 = err12;
e2 = err6;
h = 1
T = size(e1,1);
% Define the loss differential
d = e1.^2 - e2.^2;
% Ralculate the variance of the loss differential, taking into account
% autocorrelation.
dMean = mean(d);
gamma0 = var(d);
if h > 1
    gamma = zeros(h-1,1);
    for i = 1:h-1
        gamma(i) = ( d(1+i:T)' * d(1:T-i) ) ./ T; % bugfix
    end
    varD = gamma0 + 2*sum(gamma);
else
    varD = gamma0;
end
% Retrieve the diebold mariano statistic DM ~N(0,1)
DM = dMean / sqrt ( (1/T)*varD );  

%P_VALUE is calculated    
    p_value = 2*tcdf(-abs(DM),T-1)
    
    
%% Plot CSSED

serr1 = err_RWnodrift.^2;
serr2 = err_RWdriftrolling.^2;
serr3 = err_RWdriftexpanding.^2;
serr4 = err_AR1rolling.^2;
serr5 = err_AR1expanding.^2;
serr6 = err_Linearrolling.^2;
serr7 = err_Linearexpanding.^2;
serr8 = err_ARX1rolling.^2;
serr9 = err_ARX1expanding.^2;
serr10 = err_LSTMbetasrolling'.^2;
serr11 = err_LSTMbetasexpanding'.^2;
serr12 = err_MCrolling.^2;
serr13 = err_MCexpanding.^2;
serr14 = err_LSTMcurretrolling'.^2;
serr15 = err_LSTMcurretexpanding'.^2;
serr16 = err_LSTMbigdatarolling'.^2;
serr17 = err_LSTMbigdataexpanding'.^2;


% første ledd RW
CSSED_1 = cumsum(serr1 - serr2);

CSSED_2 = cumsum(serr1 - serr3);

CSSED_3 = cumsum(serr1 - serr4);

CSSED_4 = cumsum(serr1 - serr5);

CSSED_5 = cumsum(serr1 - serr6);

CSSED_6 = cumsum(serr1 - serr7);

CSSED_7 = cumsum(serr1 - serr8);

CSSED_8 = cumsum(serr1 - serr9);

CSSED_9 = cumsum(serr1 - serr10);

CSSED_10 = cumsum(serr1 - serr11);

CSSED_11 = cumsum(serr1 - serr12);

CSSED_12 = cumsum(serr1 - serr13);

CSSED_13 = cumsum(serr1 - serr14);

CSSED_14 = cumsum(serr1 - serr15);

CSSED_15 = cumsum(serr1 - serr16);

CSSED_16 = cumsum(serr1 - serr17);

figure
set(gcf,'color','w');

subplot(4,1,1)
plot(CSSED_1, 'LineWidth', 1.1)
hold on
plot(CSSED_2,'-.', 'LineWidth', 1.1)
hold on
yline(0)
hold off
title("a) Random Walk with drift")
legend('Rolling Window','Expanding Window','Orientation','horizontal')
ylabel("CSSED")
xlabel("Year")
ytickangle(90)
xticks(11:24:182)
xticklabels([2007:2:2021])
xlim([1 182])
box off

subplot(4,1,2)
plot(CSSED_3, 'LineWidth', 1.1)
hold on
plot(CSSED_4,'-.', 'LineWidth', 1.1)
hold on
yline(0)
hold off
title({'';'';"b) AR(1)"})
legend('Rolling Window','Expanding Window','Orientation','horizontal')
ylabel("CSSED")
xlabel("Year")
ytickangle(90)
xticks(11:24:182)
xticklabels([2007:2:2021])
xlim([1 182])
box off

subplot(4,1,3)
plot(CSSED_5, 'LineWidth', 1.1)
hold on
plot(CSSED_6,'-.', 'LineWidth', 1.1)
hold on
yline(0)
hold off
title({'';'';"c) Linear Regression"})
legend('Rolling Window','Expanding Window','Orientation','horizontal')
ylabel("CSSED")
xlabel("Year")
ytickangle(90)
xticks(11:24:182)
xticklabels([2007:2:2021])
xlim([1 182])
box off

subplot(4,1,4)
plot(CSSED_7, 'LineWidth', 1.1)
hold on
plot(CSSED_8,'-.', 'LineWidth', 1.1)
hold on
yline(0)
hold off
title({'';'';"d) ARX(1)"})
legend('Rolling Window','Expanding Window','Orientation','horizontal')
ylabel("CSSED")
xlabel("Year")
ytickangle(90)
xticks(11:24:182)
xticklabels([2007:2:2021])
xlim([1 182])
box off


figure
set(gcf,'color','w');

subplot(4,1,1)
plot(CSSED_9, 'LineWidth', 1.1)
hold on
plot(CSSED_10,'-.', 'LineWidth', 1.1)
hold on
yline(0)
hold off
title("e) LSTM Network with Level, Slope and Curvature")
legend('Rolling Window','Expanding Window','Orientation','horizontal')
ylabel("CSSED")
xlabel("Year")
ytickangle(90)
xticks(11:24:182)
xticklabels([2007:2:2021])
xlim([1 182])
box off

subplot(4,1,2)
plot(CSSED_11, 'LineWidth', 1.1)
hold on
plot(CSSED_12,'-.', 'LineWidth', 1.1)
hold on
yline(0)
hold off
title({'';'';"f) Model Combination"})
legend('Rolling Window','Expanding Window','Orientation','horizontal')
ylabel("CSSED")
xlabel("Year")
ytickangle(90)
xticks(11:24:182)
xticklabels([2007:2:2021])
xlim([1 182])
box off

subplot(4,1,3)
plot(CSSED_13, 'LineWidth', 1.1)
hold on
plot(CSSED_14,'-.', 'LineWidth', 1.1)
hold on
yline(0)
hold off
title({'';'';"g) LSTM Network with Currency Return"})
legend('Rolling Window','Expanding Window','Orientation','horizontal')
ylabel("CSSED")
xlabel("Year")
ytickangle(90)
xticks(11:24:182)
xticklabels([2007:2:2021])
xlim([1 182])
box off

subplot(4,1,4)
plot(CSSED_15, 'LineWidth', 1.1)
hold on
plot(CSSED_16,'-.', 'LineWidth', 1.1)
hold on
yline(0)
hold off
title({'';'';"h) LSTM Network with big dataset"})
legend('Rolling Window','Expanding Window','Orientation','horizontal')
ylabel("CSSED")
xlabel("Year")
ytickangle(90)
xticks(11:24:182)
xticklabels([2007:2:2021])
xlim([1 182])
box off



%% Revenue graph
actCurRet = CurrencyReturn(61:end);
wealth1 = zeros(183,1);
wealth1(1) = 1;
yhat1 = yhat_RWnodrift;
for i = 2:183
    if yhat1(i-1) > 0 && actCurRet(i-1) > 0 
        wealth1(i) = wealth1(i-1)*(1 + actCurRet(i-1));
    end
    if yhat1(i-1) > 0 && actCurRet(i-1) < 0
        wealth1(i) = wealth1(i-1)*(1 + actCurRet(i-1));
    end
    if yhat1(i-1) < 0 && actCurRet(i-1) < 0
        wealth1(i) = wealth1(i-1)*(1 + abs(actCurRet(i-1)));
    end
    if yhat1(i-1) < 0 && actCurRet(i-1) > 0
        wealth1(i) = wealth1(i-1)*(1 + (-1)*actCurRet(i-1));
    end
    if yhat1(i-1) == 0
        wealth1(i) = wealth1(i-1);
    end
end

plot(wealth1);

% Revenue graph
actCurRet = CurrencyReturn(61:end);
wealth2 = zeros(183,1);
wealth2(1) = 1;
yhat2 = yhat_RWdriftrolling;
for i = 2:183
    if yhat2(i-1) > 0 && actCurRet(i-1) > 0 
        wealth2(i) = wealth2(i-1)*(1 + actCurRet(i-1));
    end
    if yhat2(i-1) > 0 && actCurRet(i-1) < 0
        wealth2(i) = wealth2(i-1)*(1 + actCurRet(i-1));
    end
    if yhat2(i-1) < 0 && actCurRet(i-1) < 0
        wealth2(i) = wealth2(i-1)*(1 + abs(actCurRet(i-1)));
    end
    if yhat2(i-1) < 0 && actCurRet(i-1) > 0
        wealth2(i) = wealth2(i-1)*(1 + (-1)*actCurRet(i-1));
    end
    if yhat2(i-1) == 0
        wealth2(i) = wealth2(i-1);
    end
end

plot(wealth2);


% Revenue graph
actCurRet = CurrencyReturn(61:end);
wealth3 = zeros(183,1);
wealth3(1) = 1;
yhat3 = yhat_RWdriftexpanding;
for i = 2:183
    if yhat3(i-1) > 0 && actCurRet(i-1) > 0 
        wealth3(i) = wealth3(i-1)*(1 + actCurRet(i-1));
    end
    if yhat3(i-1) > 0 && actCurRet(i-1) < 0
        wealth3(i) = wealth3(i-1)*(1 + actCurRet(i-1));
    end
    if yhat3(i-1) < 0 && actCurRet(i-1) < 0
        wealth3(i) = wealth3(i-1)*(1 + abs(actCurRet(i-1)));
    end
    if yhat3(i-1) < 0 && actCurRet(i-1) > 0
        wealth3(i) = wealth3(i-1)*(1 + (-1)*actCurRet(i-1));
    end
    if yhat3(i-1) == 0
        wealth3(i) = wealth3(i-1);
    end
end

plot(wealth3);


% Revenue graph
actCurRet = CurrencyReturn(61:end);
wealth4 = zeros(183,1);
wealth4(1) = 1;
yhat4 = yhat_AR1rolling;
for i = 2:183
    if yhat4(i-1) > 0 && actCurRet(i-1) > 0 
        wealth4(i) = wealth4(i-1)*(1 + actCurRet(i-1));
    end
    if yhat4(i-1) > 0 && actCurRet(i-1) < 0
        wealth4(i) = wealth4(i-1)*(1 + actCurRet(i-1));
    end
    if yhat4(i-1) < 0 && actCurRet(i-1) < 0
        wealth4(i) = wealth4(i-1)*(1 + abs(actCurRet(i-1)));
    end
    if yhat4(i-1) < 0 && actCurRet(i-1) > 0
        wealth4(i) = wealth4(i-1)*(1 + (-1)*actCurRet(i-1));
    end
    if yhat4(i-1) == 0
        wealth4(i) = wealth4(i-1);
    end
end

plot(wealth4);


% Revenue graph
actCurRet = CurrencyReturn(61:end);
wealth5 = zeros(183,1);
wealth5(1) = 1;
yhat5 = yhat_AR1expanding;
for i = 2:183
    if yhat5(i-1) > 0 && actCurRet(i-1) > 0 
        wealth5(i) = wealth5(i-1)*(1 + actCurRet(i-1));
    end
    if yhat5(i-1) > 0 && actCurRet(i-1) < 0
        wealth5(i) = wealth5(i-1)*(1 + actCurRet(i-1));
    end
    if yhat5(i-1) < 0 && actCurRet(i-1) < 0
        wealth5(i) = wealth5(i-1)*(1 + abs(actCurRet(i-1)));
    end
    if yhat5(i-1) < 0 && actCurRet(i-1) > 0
        wealth5(i) = wealth5(i-1)*(1 + (-1)*actCurRet(i-1));
    end
    if yhat5(i-1) == 0
        wealth5(i) = wealth5(i-1);
    end
end

plot(wealth5);


% Revenue graph
actCurRet = CurrencyReturn(61:end);
wealth6 = zeros(183,1);
wealth6(1) = 1;
yhat6 = yhat_Linearrolling;
for i = 2:183
    if yhat6(i-1) > 0 && actCurRet(i-1) > 0 
        wealth6(i) = wealth6(i-1)*(1 + actCurRet(i-1));
    end
    if yhat6(i-1) > 0 && actCurRet(i-1) < 0
        wealth6(i) = wealth6(i-1)*(1 + actCurRet(i-1));
    end
    if yhat6(i-1) < 0 && actCurRet(i-1) < 0
        wealth6(i) = wealth6(i-1)*(1 + abs(actCurRet(i-1)));
    end
    if yhat6(i-1) < 0 && actCurRet(i-1) > 0
        wealth6(i) = wealth6(i-1)*(1 + (-1)*actCurRet(i-1));
    end
    if yhat6(i-1) == 0
        wealth6(i) = wealth6(i-1);
    end
end

plot(wealth6);


% Revenue graph
actCurRet = CurrencyReturn(61:end);
wealth7 = zeros(183,1);
wealth7(1) = 1;
yhat7 = yhat_Linearexpanding;
for i = 2:183
    if yhat7(i-1) > 0 && actCurRet(i-1) > 0 
        wealth7(i) = wealth7(i-1)*(1 + actCurRet(i-1));
    end
    if yhat7(i-1) > 0 && actCurRet(i-1) < 0
        wealth7(i) = wealth7(i-1)*(1 + actCurRet(i-1));
    end
    if yhat7(i-1) < 0 && actCurRet(i-1) < 0
        wealth7(i) = wealth7(i-1)*(1 + abs(actCurRet(i-1)));
    end
    if yhat7(i-1) < 0 && actCurRet(i-1) > 0
        wealth7(i) = wealth7(i-1)*(1 + (-1)*actCurRet(i-1));
    end
    if yhat7(i-1) == 0
        wealth7(i) = wealth7(i-1);
    end
end

plot(wealth7);


% Revenue graph
actCurRet = CurrencyReturn(61:end);
wealth8 = zeros(183,1);
wealth8(1) = 1;
yhat8 = yhat_ARX1rolling;
for i = 2:183
    if yhat8(i-1) > 0 && actCurRet(i-1) > 0 
        wealth8(i) = wealth8(i-1)*(1 + actCurRet(i-1));
    end
    if yhat8(i-1) > 0 && actCurRet(i-1) < 0
        wealth8(i) = wealth8(i-1)*(1 + actCurRet(i-1));
    end
    if yhat8(i-1) < 0 && actCurRet(i-1) < 0
        wealth8(i) = wealth8(i-1)*(1 + abs(actCurRet(i-1)));
    end
    if yhat8(i-1) < 0 && actCurRet(i-1) > 0
        wealth8(i) = wealth8(i-1)*(1 + (-1)*actCurRet(i-1));
    end
    if yhat8(i-1) == 0
        wealth8(i) = wealth8(i-1);
    end
end

plot(wealth8);


% Revenue graph
actCurRet = CurrencyReturn(61:end);
wealth9 = zeros(183,1);
wealth9(1) = 1;
yhat9 = yhat_ARX1expanding;
for i = 2:183
    if yhat9(i-1) > 0 && actCurRet(i-1) > 0 
        wealth9(i) = wealth9(i-1)*(1 + actCurRet(i-1));
    end
    if yhat9(i-1) > 0 && actCurRet(i-1) < 0
        wealth9(i) = wealth9(i-1)*(1 + actCurRet(i-1));
    end
    if yhat9(i-1) < 0 && actCurRet(i-1) < 0
        wealth9(i) = wealth9(i-1)*(1 + abs(actCurRet(i-1)));
    end
    if yhat9(i-1) < 0 && actCurRet(i-1) > 0
        wealth9(i) = wealth9(i-1)*(1 + (-1)*actCurRet(i-1));
    end
    if yhat9(i-1) == 0
        wealth9(i) = wealth9(i-1);
    end
end

plot(wealth9);


% Revenue graph
actCurRet = CurrencyReturn(61:end);
wealth10 = zeros(183,1);
wealth10(1) = 1;
yhat10 = yhat_LSTMbetasrolling';
for i = 2:183
    if yhat10(i-1) > 0 && actCurRet(i-1) > 0 
        wealth10(i) = wealth10(i-1)*(1 + actCurRet(i-1));
    end
    if yhat10(i-1) > 0 && actCurRet(i-1) < 0
        wealth10(i) = wealth10(i-1)*(1 + actCurRet(i-1));
    end
    if yhat10(i-1) < 0 && actCurRet(i-1) < 0
        wealth10(i) = wealth10(i-1)*(1 + abs(actCurRet(i-1)));
    end
    if yhat10(i-1) < 0 && actCurRet(i-1) > 0
        wealth10(i) = wealth10(i-1)*(1 + (-1)*actCurRet(i-1));
    end
    if yhat10(i-1) == 0
        wealth10(i) = wealth10(i-1);
    end
end

plot(wealth10);


% Revenue graph
actCurRet = CurrencyReturn(61:end);
wealth11 = zeros(183,1);
wealth11(1) = 1;
yhat11 = yhat_LSTMbetasexpanding';
for i = 2:183
    if yhat11(i-1) > 0 && actCurRet(i-1) > 0 
        wealth11(i) = wealth11(i-1)*(1 + actCurRet(i-1));
    end
    if yhat11(i-1) > 0 && actCurRet(i-1) < 0
        wealth11(i) = wealth11(i-1)*(1 + actCurRet(i-1));
    end
    if yhat11(i-1) < 0 && actCurRet(i-1) < 0
        wealth11(i) = wealth11(i-1)*(1 + abs(actCurRet(i-1)));
    end
    if yhat11(i-1) < 0 && actCurRet(i-1) > 0
        wealth11(i) = wealth11(i-1)*(1 + (-1)*actCurRet(i-1));
    end
    if yhat11(i-1) == 0
        wealth11(i) = wealth11(i-1);
    end
end

plot(wealth11);


% Revenue graph
actCurRet = CurrencyReturn(61:end);
wealth12 = zeros(183,1);
wealth12(1) = 1;
yhat12 = yhat_MCrolling;
for i = 2:183
    if yhat12(i-1) > 0 && actCurRet(i-1) > 0 
        wealth12(i) = wealth12(i-1)*(1 + actCurRet(i-1));
    end
    if yhat12(i-1) > 0 && actCurRet(i-1) < 0
        wealth12(i) = wealth12(i-1)*(1 + actCurRet(i-1));
    end
    if yhat12(i-1) < 0 && actCurRet(i-1) < 0
        wealth12(i) = wealth12(i-1)*(1 + abs(actCurRet(i-1)));
    end
    if yhat12(i-1) < 0 && actCurRet(i-1) > 0
        wealth12(i) = wealth12(i-1)*(1 + (-1)*actCurRet(i-1));
    end
    if yhat12(i-1) == 0
        wealth12(i) = wealth12(i-1);
    end
end

plot(wealth12);


% Revenue graph
actCurRet = CurrencyReturn(61:end);
wealth13 = zeros(183,1);
wealth13(1) = 1;
yhat13 = yhat_MCexpanding;
for i = 2:183
    if yhat13(i-1) > 0 && actCurRet(i-1) > 0 
        wealth13(i) = wealth13(i-1)*(1 + actCurRet(i-1));
    end
    if yhat13(i-1) > 0 && actCurRet(i-1) < 0
        wealth13(i) = wealth13(i-1)*(1 + actCurRet(i-1));
    end
    if yhat13(i-1) < 0 && actCurRet(i-1) < 0
        wealth13(i) = wealth13(i-1)*(1 + abs(actCurRet(i-1)));
    end
    if yhat13(i-1) < 0 && actCurRet(i-1) > 0
        wealth13(i) = wealth13(i-1)*(1 + (-1)*actCurRet(i-1));
    end
    if yhat13(i-1) == 0
        wealth13(i) = wealth13(i-1);
    end
end

plot(wealth13);


% Revenue graph
actCurRet = CurrencyReturn(61:end);
wealth14 = zeros(183,1);
wealth14(1) = 1;
yhat14 = yhat_LSTMcurretrolling';
for i = 2:183
    if yhat14(i-1) > 0 && actCurRet(i-1) > 0 
        wealth14(i) = wealth14(i-1)*(1 + actCurRet(i-1));
    end
    if yhat14(i-1) > 0 && actCurRet(i-1) < 0
        wealth14(i) = wealth14(i-1)*(1 + actCurRet(i-1));
    end
    if yhat14(i-1) < 0 && actCurRet(i-1) < 0
        wealth14(i) = wealth14(i-1)*(1 + abs(actCurRet(i-1)));
    end
    if yhat14(i-1) < 0 && actCurRet(i-1) > 0
        wealth14(i) = wealth14(i-1)*(1 + (-1)*actCurRet(i-1));
    end
    if yhat14(i-1) == 0
        wealth14(i) = wealth14(i-1);
    end
end

plot(wealth14);


% Revenue graph
actCurRet = CurrencyReturn(61:end);
wealth15 = zeros(183,1);
wealth15(1) = 1;
yhat15 = yhat_LSTMcurretexpanding';
for i = 2:183
    if yhat15(i-1) > 0 && actCurRet(i-1) > 0 
        wealth15(i) = wealth15(i-1)*(1 + actCurRet(i-1));
    end
    if yhat15(i-1) > 0 && actCurRet(i-1) < 0
        wealth15(i) = wealth15(i-1)*(1 + actCurRet(i-1));
    end
    if yhat15(i-1) < 0 && actCurRet(i-1) < 0
        wealth15(i) = wealth15(i-1)*(1 + abs(actCurRet(i-1)));
    end
    if yhat15(i-1) < 0 && actCurRet(i-1) > 0
        wealth15(i) = wealth15(i-1)*(1 + (-1)*actCurRet(i-1));
    end
    if yhat15(i-1) == 0
        wealth15(i) = wealth15(i-1);
    end
end

plot(wealth15);


% Revenue graph
actCurRet = CurrencyReturn(61:end);
wealth16 = zeros(183,1);
wealth16(1) = 1;
yhat16 = yhat_LSTMbigdatarolling';
for i = 2:183
    if yhat16(i-1) > 0 && actCurRet(i-1) > 0 
        wealth16(i) = wealth16(i-1)*(1 + actCurRet(i-1));
    end
    if yhat16(i-1) > 0 && actCurRet(i-1) < 0
        wealth16(i) = wealth16(i-1)*(1 + actCurRet(i-1));
    end
    if yhat16(i-1) < 0 && actCurRet(i-1) < 0
        wealth16(i) = wealth16(i-1)*(1 + abs(actCurRet(i-1)));
    end
    if yhat16(i-1) < 0 && actCurRet(i-1) > 0
        wealth16(i) = wealth16(i-1)*(1 + (-1)*actCurRet(i-1));
    end
    if yhat16(i-1) == 0
        wealth16(i) = wealth16(i-1);
    end
end

plot(wealth16);


% Revenue graph
actCurRet = CurrencyReturn(61:end);
wealth17 = zeros(183,1);
wealth17(1) = 1;
yhat17 = yhat_LSTMbigdataexpanding';
for i = 2:183
    if yhat17(i-1) > 0 && actCurRet(i-1) > 0 
        wealth17(i) = wealth17(i-1)*(1 + actCurRet(i-1));
    end
    if yhat17(i-1) > 0 && actCurRet(i-1) < 0
        wealth17(i) = wealth17(i-1)*(1 + actCurRet(i-1));
    end
    if yhat17(i-1) < 0 && actCurRet(i-1) < 0
        wealth17(i) = wealth17(i-1)*(1 + abs(actCurRet(i-1)));
    end
    if yhat17(i-1) < 0 && actCurRet(i-1) > 0
        wealth17(i) = wealth17(i-1)*(1 + (-1)*actCurRet(i-1));
    end
    if yhat17(i-1) == 0
        wealth17(i) = wealth17(i-1);
    end
end

plot(wealth17);

%% Plot of all wealth graphs

figure
set(gcf,'color','w');
plot(wealth1, 'k', 'LineWidth', 1.1)
hold on
plot(wealth2, 'color', '#D95319', 'LineWidth', 1.1)
hold on
plot(wealth3,'-.', 'color', '#D95319', 'LineWidth', 1.1)
hold on
plot(wealth4, 'color', '#77AC30', 'LineWidth', 1.1)
hold on
plot(wealth5,'-.', 'color', '#77AC30', 'LineWidth', 1.1)
hold on
plot(wealth6, 'm', 'LineWidth', 1.1)
hold on
plot(wealth7,'-.', 'color', 'm', 'LineWidth', 1.1)
hold on
plot(wealth8, 'y', 'LineWidth', 1.1)
hold on
plot(wealth9, '-.', 'color', 'y', 'LineWidth', 1.1)
hold on
plot(wealth10, 'c', 'LineWidth', 1.1)
hold on
plot(wealth11,'-.', 'color', 'c', 'LineWidth', 1.1)
hold on
plot(wealth12, 'g', 'LineWidth', 1.1)
hold on
plot(wealth13,'-.', 'color', 'g', 'LineWidth', 1.1)
hold on
plot(wealth14, 'r', 'LineWidth', 1.1)
hold on
plot(wealth15,'-.', 'color', 'r', 'LineWidth', 1.1)
hold on
plot(wealth16, 'b', 'LineWidth', 1.1)
hold on
plot(wealth17,'-.', 'color', 'b', 'LineWidth', 1.1)
hold off
legend('RW without drift','RW drift rolling','RW drift expanding','AR(1) rolling', 'AR(1) expanding', 'Linear rolling', 'Linear expanding', 'ARX(1) rolling', 'ARX(1) expanding', 'AI betas rolling', 'AI betas expanding', 'MC rolling', 'MC expanding', 'AI currency return rolling', 'AI currency return expanding', 'AI big dataset rolling', 'AI big dataset expanding')
ylabel("Wealth")
xlabel("Year")
ytickangle(90)
xticks(11:24:182)
xticklabels([2007:2:2021])
xlim([1 182])
box off




