%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Empirical Financial Econometrics Matlab Class Code part 3
% Updated Spring 2021
% This Matlab code facilitates to ASSN 5-6. 
% Copyright at Chaoyi Chen (BCE & MNB)
% All rights reserved. For use by registered students only.
% Please do NOT distribute without express written consent.
% Email: chaoyi.chen@uni-corvinus.hu
% Web: www.chenchaoyi.com
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
%{
Folloiwng lines of codes are from previous ASSNs as we need part of their results first
%}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
%%
clear all 
warning('off') 

%%
% Load data
load YMOPD1901 

year = YMOPD(:,1); 
price = YMOPD(:,2);
lagprice = lagmatrix(price,1);
n=size(year,1);

% Define log versions 
p = log(price);
p = p(2:end);
lp = log(lagprice);
lp = lp(2:end);

% Define ratios and returns
y = p - lp; 

pMax = 1; 
qMax = 3; 
LogL = zeros(pMax, qMax);
SumPQ = LogL;
for p = 1:pMax
    for q = 1:qMax
        ToEstMdl = arima(p,0,q); % construct ARMA(p,q) model
        [Est, ~, LogL(p,q)] = estimate(ToEstMdl, y, 'Display','off'); % estimate ARMA(p,q) model
        SumPQ(p,q) = p+q;
    end
end

logL = reshape(LogL, pMax*qMax,1); % reshape our log-likelihood for each pair of (p,q)
numParams = reshape(SumPQ,pMax*qMax,1)+1; % reshape number of parameters used in each try
bic = aicbic(logL,numParams,size(y,1)); % compute BIC for each pair
BIC = reshape(bic, pMax, qMax); % transform back to pMax by qMax matrix to help us locate the optimal (p,q)
[bestP,bestQ] = find(BIC == min(bic)); % find optimal p and q

ToEstMdl_opt = arima(bestP,0,bestQ);
[Est_opt, ~, info_opt] = estimate(ToEstMdl_opt, y, 'Display','off');

ARMA_resid = infer(Est_opt,y); % collect in sample residual. 
ARMA_resid_squared = ARMA_resid.^2;
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
%{
Now begin with the new study
%}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
%%
%%%%%%%%%%%%%%%%%%%%%%%%% GARCH model selection %%%%%%%%%%%%%%%%%%%%%%%%%
% As we find evidence of ARCH effect, we proceed to use GARCH model modling
% the conditional volatility. 

% Similar to ARMA, set largest lags for both ARCH and GARCH terms and estimate all models.
pMax_garch=3;
qMax_garch=3;
LogL_garch=zeros(pMax_garch,qMax_garch);
SumPQ_garch=LogL_garch;
for p= 1:pMax_garch
for q = 1:qMax_garch
ToEstMdl_garch=garch(p,q);
[Est_garch,estco_garch,LogL_garch(p,q)]=estimate(ToEstMdl_garch,ARMA_resid,'Display','off');
SumPQ_garch(p,q)=p+q;
end
end

% select the optimal pq by BIC
logL_garch=reshape(LogL_garch,pMax_garch*qMax_garch,1);
numParams_garch=reshape(SumPQ_garch,pMax_garch*qMax_garch,1)+1;
bic_garch=aicbic(logL_garch,numParams_garch,size(ARMA_resid,1));
BIC_garch=reshape(bic_garch,pMax_garch,qMax_garch);
[bestP_garch,bestQ_garch]=find(BIC_garch==min(bic_garch));

% after we find the optimal (p,q), we proceed to construct the optimal GARCH model and
% estimate the model
ToEstMdl_opt_garch = garch(bestP_garch,bestQ_garch);
[Est_opt_garch, ~, info_opt_garch] = estimate(ToEstMdl_opt_garch, ARMA_resid, 'Display','off');

% Print results of in sample results of optimal GARCH mdoel
%---------TABLES-----------%
diary results_GARCH_is; 
    disp('---------------------------');
    disp('GARCH in sample results');
    disp('---------------------------');
    disp(Est_opt_garch);
    disp('---------------------------')
diary off;

%%
%%%%%%%%%%%%%%%%%%%%%%%%% In sample analysis %%%%%%%%%%%%%%%%%%%%%%%%%

% get the conditional varianc of ARMA residual
ARMA_resid_v = infer(Est_opt_garch,ARMA_resid);

% Plot the infered conditional variance by GARCH versus the realized the variance
set(0,'DefaulttextInterpreter','latex');
plot(year(2:end),ARMA_resid_v,  'LineWidth',2);
hold on
plot(year(2:end),ARMA_resid_squared, 'Color',[.7,.7,.7], 'LineWidth',2);
hold on
legend({'GARCH', 'Realized'},'Interpreter','latex','location','NW','FontSize',15);
xlabel('Year','FontSize',10);
ylabel('Conditional Variance','FontSize',10);
title('Infered Conditional Variance by Optimal GARCH versus the Realized the Variance','FontSize',15)
grid on;

%%
%%%%%%%%%%%%%%%%%%%%%%%%% Out-of-sample analysis %%%%%%%%%%%%%%%%%%%%%%%%%
% In the following, we aim to evaluate the optimal GARCH model's OOS preformance

% Similar to OOS ARMA, we first divide the sample into two part, training and testing
training_garch = 40; % set initial training period for recursive or rolling window
estMethod_garch = 'recursivewindow'; % specify OOS method. Option: 'recursivewindow' or 'rollingwindow'

fitted_y_garch = zeros(size(ARMA_resid,1)-training_garch,1); % store fitted conditional var of GARCH

for t=1:size(ARMA_resid,1)-training_garch
    
    % Specify training data
    if strcmpi(estMethod_garch,'rollingwindow')
       Esty_garch = ARMA_resid(t:training_garch+t-1);
    else
       if strcmpi(estMethod_garch,'recursivewindow') 
           Esty_garch = ARMA_resid(1:training_garch+t-1);
       else
           warning('either rollingwindow or recursivewindow ha to be choosen')
       end
    end
    
    % Specify testing data
    Outy_garch = ARMA_resid(training_garch+t);
    
    % Estimate optimal ARMA model
    [Est_opt_oos_garch, ~, info_opt_oos_garch] = estimate(ToEstMdl_opt_garch, Esty_garch, 'Display','off');
    
    % Forecast next period return from ARMA model
    fitted_y_garch(t) = forecast(Est_opt_oos_garch,1,'Y0',Esty_garch);
end    

% Plot the forecast
figure
h1 = plot(ARMA_resid_v,'Color',[.7,.7,.7],'LineWidth',2);
hold on
h2 = plot(training_garch+1:size(ARMA_resid_v,1),fitted_y_garch,'b','LineWidth',2);
legend([h1 h2],'Realized','Forecast by GARCH','Location','NorthWest');
title('Out-of-sample Forecasts for Conditional Variance');
grid on
hold off

%%
%%%%%%%%%%%%%%%%%%%%%%%%% Dynamic mean-variance portfolio analysis %%%%%%%%%%%%%%%%%%%%%%%%%
%{
 In this study, we want to build a dynamic mean-variance portfolio. To
 simplify the study, we impose some regularity assumptions on the assets.

 1). We assume risk free asset has the same return (R_f = 0.01) over years
 2). We also assume our investor is risk averse with k=1 (see lecture note session 1 page 47 for more reference)
 3). we do NOT allow for short risky asset (that is w_{t} has to
 be nonnegative) and leverage on ivestment (e.g. w_{t}>1) so w_{t} has to be between 0 to 1. However, feel free to include the short operation at
 home.

 We will compare our dynamic portfolio with w = 0, 1/2, 1 by the out-of-sample method
 Note that, our dynamic weight will vary each period as optimal w_{t}=(1/k)[(E_t(R_{t+1})-R_f)/var_t(R_{t+1})] (see lecture note session 1 page 49 for more reference)

 As the investor is risk averse, it is meaningless to only look at the
 return. We will evaluate the utility gained in each portfolio as well by
 using the mean-variance utility function we defined in lecture note (session 1) page 48
%}

% We use optimal ARMA [(bestP,bestQ)] to forecast R_{t+1} and use GARCG [(bestP_garch,bestQ_garch)] to
% foreacast var_t(R_{t+1})
% ToEstMdl_opt_garch = garch(1,0); 

riskfree = 0.03; % assume risk free asset has the same return = 0.01
rf_vector = riskfree*ones(size(y,1),1);
k = 0.5; % assume investor is risk averse with k
training_portfolio = 60; % unlike to compare static portfolio (we can compare all sample periods), we have to lose some observations to compute the dynamic weight and c=only compare testing period
estMethod_portfolio = 'recursivewindow'; % method to select dynamic weight

return_0 = zeros(size(y,1)-training_portfolio,1); % store static fortfolio return with w = 0 for all t
return_half = zeros(size(y,1)-training_portfolio,1); % store static fortfolio return with w = 1/2 for all t
return_1 = zeros(size(y,1)-training_portfolio,1); % store static fortfolio return with w = 1 for all t

return_0_forecast = zeros(size(y,1)-training_portfolio,1); % store static fortfolio forecast with w = 0 for all t
return_half_forecast = zeros(size(y,1)-training_portfolio,1); % store static fortfolio forecast with w = 1/2 for all t
return_1_forecast = zeros(size(y,1)-training_portfolio,1); % store static fortfolio forecast with w = 1 for all t
return_dynamic_forecast = zeros(size(y,1)-training_portfolio,1); % store static fortfolio forecast with w = 1 for all t

fitted_y_portfolio = zeros(size(y,1)-training_portfolio,1); % store ARMA estimate of E_t(R_{t+1})
fitted_y_garch_portfolio = zeros(size(y,1)-training_portfolio,1); % store GARCH estimate of Var_t(R_{t+1})

w_optimal = zeros(size(y,1)-training_portfolio,1); % store dynamic optimal weight for each period
return_dynamic = zeros(size(y,1)-training_portfolio,1); % store dynamic fortfolio return with w = w_optimal(t)

u_0 = zeros(size(y,1)-training_portfolio,1); % store utility gained from static fortfolio with w = 0 for all t
u_half = zeros(size(y,1)-training_portfolio,1); % store utility gained from static fortfolio with w = 1/2 for all t
u_1 = zeros(size(y,1)-training_portfolio,1); % store utility gained from static fortfolio with w = 1 for all t
u_dynamic = zeros(size(y,1)-training_portfolio,1); % store utility gained from static fortfolio with dynamic w for all t

for t=1:size(y,1)-training_portfolio
    
    % Specify training data
    if strcmpi(estMethod_portfolio,'rollingwindow')
       Esty_portfolio = y(t:training_portfolio+t-1);
    else
       if strcmpi(estMethod_portfolio,'recursivewindow') 
           Esty_portfolio = y(1:training_portfolio+t-1);
       else
           warning('either rollingwindow or recursivewindow ha to be choosen')
       end
    end
    
    % Specify testing data/ next period true return/ vary on risky asset
    Outy_portfolio = y(training_portfolio+t);
    
    %%%%%%%%%%%%%%%%%%%%%%%% Static (weight) portfolio return %%%%%%%%%%%%%%%%%%%%%%%%
    % forecast next period return for each portfolio
    return_0_forecast(t) = riskfree; 
    return_half_forecast(t) = 0.5*mean(Esty_portfolio)+0.5*riskfree;
    return_1_forecast(t) = mean(Esty_portfolio); 
    
    % realized next period return for each portfolio
    return_0(t) = riskfree; 
    return_half(t) = 0.5*Outy_portfolio+0.5*riskfree;
    return_1(t) = Outy_portfolio(end); 
    
    %%%%%%%%%%%%%%%%%%%%%%%% Dynamic portfolio return %%%%%%%%%%%%%%%%%%%%%%%%
    %{
     We solve the dynamic portfolio in the following steps:
     1). We use ARMA estimate E_t(R_{t+1}) and collect the training period
     residual
     2). We use GARCH estimate var_t(R_{t+1}). 
     3). Then, apply w_{t}=(1/k)[(E_t(R_{t+1})-R_f)/var_t(R^{p}_{t+1})] to comput
     each period weight
     4). Lastly, we obtain the dynamic portfolio return_{t+1} =
     w_{t}*R_{t+1}+(1-w_{t})*R_f
    %}
    
    % Dynamic weight selection/ obtain value of weighting function numerator
    [Est_opt_oos_portfolio, ~, ~] = estimate(ToEstMdl_opt, Esty_portfolio, 'Display','off');
    fitted_y_portfolio(t) = forecast(Est_opt_oos_portfolio,1,'Y0',Esty_portfolio); % E_t[R_(t+1)] estimate
    ARMA_resid_portfolio = infer(Est_opt_oos_portfolio,Esty_portfolio); 
    
    % Dynamic weight selection/ btain value of weighting function denominator 
    [Est_opt_oos_garch_portfolio, ~, info_opt_oos_portfolio] = estimate(ToEstMdl_opt_garch, ARMA_resid_portfolio, 'Display','off');
    fitted_y_garch_portfolio(t) = forecast(Est_opt_oos_garch_portfolio,1,'Y0',ARMA_resid_portfolio);
    
    % Compute updated weight at period t 
    w_optimal(t) = (1/k)*(fitted_y_portfolio(t)-riskfree)/fitted_y_garch_portfolio(t);
    
    % Rule out short operation and leverage
    if w_optimal(t)<0
        w_optimal(t)=0;
    else
        if w_optimal(t)>1
        w_optimal(t)=1;
        else
           w_optimal(t)=w_optimal(t);
        end
    end
    
    % Compute forecast next period return from dynamic portfolio
    return_dynamic_forecast(t) = fitted_y_portfolio(t);
    
    % Compute realized dynamic portfolio return at t+1
    return_dynamic(t) = w_optimal(t)*Outy_portfolio+(1-w_optimal(t))*riskfree;
    
    %%%%%%%%%%%%%%%%%%%%%%%% Economic utility gained in each portfolio %%%%%%%%%%%%%%%%%%%%%%%%
    % We assume investor has utility function = return_{t+1}-k/2*var(t+1)
    % Check lecture note (session 1) page 48 for more reference
    u_0(t) = riskfree; % var=0
    u_half(t) = 0.5*riskfree+0.5*(Outy_portfolio-riskfree)-(k/2)*(1/4)*(return_half_forecast(t)-Outy_portfolio)^2;
    u_1(t) = (Outy_portfolio-riskfree)-(k/2)*(return_1_forecast(t)-Outy_portfolio)^2;
    u_dynamic(t) = (1-w_optimal(t))*riskfree+w_optimal(t)*(Outy_portfolio-riskfree)-(k/2)*(w_optimal(t))^2*(return_dynamic_forecast(t)-Outy_portfolio)^2;
end

% As we choose Training = 60, this means we initially stand at year 1961.
% At year 1961, we invest $100 on each portfolios (w=0, 1/2, 1, and dynamic)
principle = 100; % initial investment
years_invest = size(y,1)-training_portfolio;

% Portfolios start to generate profit from 1962-2015. 
% Now, we compuet the cumulative profit for each portfolio
cumu_w_0 = principle; % cumulative value of w=0 portfolio/ initial cumulative = principle
cumu_w_1 = principle; % cumulative value of w=1 portfolio/ initial cumulative = principle
cumu_w_half = principle; % cumulative value of w=1/2 portfolio/ initial cumulative = principle
cumu_w_dynamic = principle; % cumulative value of dynamic w portfolio/ initial cumulative = principle
cumu_w_0_yearly = zeros(years_invest,1); % store cumulative value for each investment year for w=0 portfolio
cumu_w_1_yearly = zeros(years_invest,1); % store cumulative value for each investment year for w=1 portfolio
cumu_w_half_yearly = zeros(years_invest,1); % store cumulative value for each investment year for w=1/2 portfolio
cumu_w_dynamic_yearly = zeros(years_invest,1); % store cumulative value for each investment year for dynamic portfolio

for i=1:years_invest
    cumu_w_0 = cumu_w_0*(1+return_0(i)); 
    cumu_w_0_yearly(i) = cumu_w_0;
    
    cumu_w_1 = cumu_w_1*(1+return_1(i)); 
    cumu_w_1_yearly(i) = cumu_w_1;
    
    cumu_w_half = cumu_w_half*(1+return_half(i)); 
    cumu_w_half_yearly(i) = cumu_w_half;
    
    cumu_w_dynamic = cumu_w_dynamic*(1+return_dynamic(i)); 
    cumu_w_dynamic_yearly(i) = cumu_w_dynamic;
end

% Plot the cumulative value series for each portfolio
set(0,'DefaulttextInterpreter','latex');
subplot(2,1,1)
plot(cumu_w_0_yearly, 'LineWidth',4);
hold on
plot(cumu_w_1_yearly, '-.', 'LineWidth',4);
hold on
plot(cumu_w_half_yearly, '-.', 'LineWidth',4);
hold on
plot(cumu_w_dynamic_yearly, '-.', 'LineWidth',4);
hold on
legend({ 'w=0', 'w=1', 'w=1/2', 'dynamic w'},'Interpreter','latex','location','NW','FontSize',15);
xlabel('Period','FontSize',10);
ylabel('Cumulative value','FontSize',10);
title('Out-of-sample Portfolio Comparison, Profit, 1962-2015','FontSize',15)
grid on;

% Print results of cumulative value of each portfolio
%---------TABLES-----------%
diary results_cumu_all; 
    disp('---------------------------');
    disp('Cumulative value of portfolio w=0');
    disp('---------------------------');
    disp(cumu_w_0);
    disp('---------------------------')
    disp('Cumulative value of portfolio w=1');
    disp('---------------------------');
    disp(cumu_w_1);
    disp('---------------------------')
    disp('Cumulative value of portfolio w=1/2');
    disp('---------------------------');
    disp(cumu_w_half);
    disp('---------------------------')
    disp('Cumulative value of dynamic weight portfolio');
    disp('---------------------------');
    disp(cumu_w_dynamic);
    disp('---------------------------')
diary off;

% Now we compute the cumulative utility gained in each testing period
cumu_u_0 = 0;
cumu_u_1 = 0;
cumu_u_half = 0;
cumu_u_dynamic = 0;
cumu_u_0_yearly = zeros(years_invest,1);
cumu_u_1_yearly = zeros(years_invest,1);
cumu_u_half_yearly = zeros(years_invest,1);
cumu_u_dynamic_yearly = zeros(years_invest,1);

for i=1:years_invest
    cumu_u_0 = cumu_u_0+u_0(i);
    cumu_u_1 = cumu_u_1+u_1(i);
    cumu_u_half = cumu_u_half+u_half(i);
    cumu_u_dynamic = cumu_u_dynamic+u_dynamic(i);
  
    cumu_u_0_yearly(i) = cumu_u_0;
    cumu_u_1_yearly(i) = cumu_u_1;
    cumu_u_half_yearly(i) = cumu_u_half;
    cumu_u_dynamic_yearly(i) = cumu_u_dynamic;
end

% Print results of cumulative utility gained from each portfolio
%---------TABLES-----------%
diary results_cumu_all_utility; 
    disp('---------------------------');
    disp('Cumulative utility gained from portfolio w=0');
    disp('---------------------------');
    disp(cumu_u_0);
    disp('---------------------------')
    disp('Cumulative utility gained from portfolio w=1');
    disp('---------------------------');
    disp(cumu_u_1);
    disp('---------------------------')
    disp('Cumulative utility gained from portfolio w=1/2');
    disp('---------------------------');
    disp(cumu_u_half);
    disp('---------------------------')
    disp('Cumulative utility gained from dynamic weight portfolio');
    disp('---------------------------');
    disp(cumu_u_dynamic);
    disp('---------------------------')
diary off;

% Plot the cumulative utilities gained in each portfolio
set(0,'DefaulttextInterpreter','latex');
subplot(2,1,2)
plot(cumu_u_0_yearly, 'LineWidth',4);
hold on
plot(cumu_u_1_yearly, '-.', 'LineWidth',4);
hold on
plot(cumu_u_half_yearly, '-.', 'LineWidth',4);
hold on
plot(cumu_u_dynamic_yearly, '-.', 'LineWidth',4);
hold on
legend({ 'w=0', 'w=1', 'w=1/2', 'dynamic w'},'Interpreter','latex','location','NW','FontSize',15);
xlabel('Period','FontSize',10);
ylabel('Cumulative utility','FontSize',10);
title('Out-of-sample Portfolio Comparison, Utility, 1962-2015','FontSize',15)
grid on;

%%%%%%%%%%%%%%%%%%%%%%%%% Home Work %%%%%%%%%%%%%%%%%%%%%%%%%
%{
1. Suppose you are a portfolio manager for a risk averse investor with k =1,
re-evaluate this study and compare both profit and utility gained from all
portfolios.
2. Suppoese you are a portfolio manager for a risk-neutral investor with k
= 0, re-evaluate this study what did you find?
3. Suppose the risk free (rf) assets is also cahnging over time, please compute
the new formular to compute the optimal weight and simulate rf using the function
rf = normrnd(0.03,0.1,size(y,1),1); Re-study the problem with your
simulated data and redefined optimal weight. 
%}
%%%%%%%%%%%%%%%%%%%%%%%%% The End %%%%%%%%%%%%%%%%%%%%%%%%%