
clear all; clc;
load data.txt;
%%Please use the spot and futures price series
spot=data(:,1);
future=data(:,2);

forecast_step=1; %%In the Management Science paper we use the forecast_step of 1. One can do with the longer forecast steps.
T=size(spot,1);
index=1;      
window=500
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Let's first calculate the residuals using VEC(2) model from OLS
%% regression for each pair of spot & future index
%% Here, we give the MRS-GARCH codes for VECM residuals. Note that in the
%% Management Science paper we do MRS-GARCH with returns directly. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    for i=1:T-window
        ret=[diff(log(spot(:,index))), diff(log(future(:,index)))]; % return of spot and future
        rolling_ret=ret(i:i+window-1,:);
        
        % Error Cointegration calculation using OLS
        log_price=[log(spot(:,index)), log(future(:,index))]; % log price of spot and future        
        rolling_log_price=log_price(i:i+window-1,:);        
        OLS_ECM=[ones(size(rolling_log_price,1),1),rolling_log_price(:,2)]\[rolling_log_price(:,1)]; %Log_s=a+b*Log_f
        ECM=rolling_log_price(:,1)-[ones(size(rolling_log_price,1),1),rolling_log_price(:,2)]*OLS_ECM;
        
        % test in Stata(R), we use VEC(2), Lag(1)_ECM for those data sets        
        y=[rolling_ret(3:end,1),rolling_ret(3:end,2)];
        x=[ones(size(y,1),1), ECM(2:end-1,:),rolling_ret(2:end-1,1),rolling_ret(2:end-1,2),...
            rolling_ret(1:end-2,1), rolling_ret(1:end-2,2)];
        coeff=x\y;
        resid{index,i}=y-x*coeff;
    end
 %%%if you want do with the return directly, using following codes
%%  for i=1:T-window
%%        ret=[100*diff(log(spot(:,index))), 100*diff(log(future(:,index)))]; % return of spot and future
%%        sample_ret{index,i}=ret(i:i+window-1,:);
%%    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% ************* Moddel 10: MRS-Scalar-BEKK-GARCH *************************
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   % for i=1:ceil((T-window)/forecast_step)
        
    %    [parameters, Ht_forecast_MRS_Scalar_BEKK, Ht_MRS_Scalar_BEKK, Rt_forecast_MRS_Scalar_BEKK, Rt_MRS_Scalar_BEKK]=...
     %       MLE_RS_Multi_Scalar_BEKK_GARCH_forecast(resid{index,(i-1)*forecast_step+1}, forecast_step);  
        
     %   Hedge_ratio_MRS_Scalar_BEKK_1(i,index)=Ht_forecast_MRS_Scalar_BEKK{1}(1,2,end)/Ht_forecast_MRS_Scalar_BEKK{1}(2,2,end);   % Hedge_ratio = cov(r1,r2)/var(r2) where r2 is the return of future index
    %end
    
    for i=1:ceil((T-window)/forecast_step)
        [parameters, Ht_forecast_MRS_DCC_GARCH, Ht_MRS_DCC_GARCH, Rt_forecast_MRS_DCC_GARCH, Rt_MRS_DCC_GARCH]=...
            MLE_RS_Multi_DCC_forecast(resid{index,(i-1)*forecast_step+1}, forecast_step);

        % form the hedge ratio and calculate the portfolio return
        Hedge_ratio_MRS_DCC_GARCH_4(i,index)=Ht_forecast_MRS_DCC_GARCH{1}(1,2,end)/Ht_forecast_MRS_DCC_GARCH{1}(2,2,end);   % Hedge_ratio = cov(r1,r2)/var(r2) where r2 is the return of future index
    end
