function forecasting_GMVPloss(n, T, spec)
%% forecasting_GMVPloss.m
% In this function, the parameters para_opt of the RLS_EF, RLS_REF or dGAS
% model are estimated; portfolio weights, portfolio return and average loss
% at the optimal parameters para_opt is predicted for a n x T return series
% r.
% The notation follows Reh et al. (2021)
% This version: 08/2021

% n = number of cross sectional units
% T = number of time units
% spec = 'RLS_EF', spec = 'RLS_REF', spec = 'dGAS'


%% load data and initialization
% load('testdata.mat');
r       = data(1:n,1:T);
y       = r(n,:);
x       = [ones(1,T); y - r(1:n-1,:)];
wgmv    = zeros(n,T);
pf_GMV  = zeros(T,1);


%% settings for optimization
optc  = optimoptions('fmincon', ...
    'GradObj', 'off', ...
    'Display','iter', ...
    'OptimalityTolerance', 1e-7, ...
    'MaxFunEvals', 10000, ...
    'MaxIter', 5000, ...
    'TolFun', 1e-16, ...
    'TolX',1e-16);

switch spec
    case 'dGAS'
        para_0 = [0, 1-1e3*eps, 0, 1-1e3*eps, 0.999];
        lb = [-inf, -1, -inf, -1, 0]; 
        ub = [inf, 1, inf, 1, 1];
    otherwise
        para_0 = [0.995, 1];
        lb = [0, 0];
        ub = [1, inf];
end


%% determination of sig2_R and rho_R
rho_0 = 0.5;
llb = -1/(n-1); 
uub = 1;
[rho_hatR] = fmincon(@ML_est_equ,rho_0,[],[],[],[],llb,uub,[],optc,r);
[~,~,SigstarR] = ML_est_equ(rho_hatR, r);
pre = 2*(1 - rho_hatR)*SigstarR(1,1);


%% estimation and evaluation
switch spec
    case 'RLS_REF'
        [para_opt] = fmincon(@rls_REF,para_0,[],[],[],[],lb,ub,[],optc, y,x,pre);
        lamhat = para_opt(1);
        gamhat = para_opt(2);
        [loss, bgmv, lo_GMV] = rls_REF(para_opt, y, x, pre);
        
    case 'RLS_EF'
        [para_opt] = fmincon(@rls_EF,para_0,[],[],[],[],lb,ub,[],optc, y,x,pre);
        lamhat = para_opt(1);
        gamhat = para_opt(2);
        [loss, bgmv, lo_GMV] = rls_EF(para_opt,y,x,pre);
        
    case 'dGAS'
        [para_opt] = fmincon(@dgas,para_0,[],[],[],[],lb,ub,[],optc, y,x,pre);
        a0hat = para_opt(3)*1e4;
        a1hat = para_opt(1)*1e4;
        b0hat = para_opt(4);
        b1hat = para_opt(2);
        kaphat = para_opt(5);
        [loss, bgmv, lo_GMV] = dgas(para_opt,y,x,pre);
end    


%% portfolio predictions
for tt = 1 : T
    
    w_          = bgmv(:,tt);
    wgmv(:,tt)  = [w_(2:n); 1-sum(w_(2:n))];
    pf_GMV(tt)  = wgmv(:,tt)'*r(:,tt);
    
end


%% saving results
fname = sprintf('%s_n%d_T%d.mat', spec, n, T);
save(fname)

end % end of main


function  [loss,b,lo]=  rls_REF(para, y, x, pre)

%% initialize and define inputs
[n, T] = size(x);
b  = zeros(n,T);
lo = zeros(T,1);

b(:,1) = 1/n*ones(n,1);
b(1,1) = 0;

lam = para(1);
gam = para(2);

En  = 0.5*eye(n-1) + 0.5*ones(n-1);
Om_ = [1 zeros(1,n-1); zeros(n-1,1), pre*En];
Om  = gam*Om_;
Om_star = Om;
Pt  = Om\eye(n);
H   = lam;
Tt  = eye(n);
In  = eye(n);


%% start recursion
for tt = 1 : T
    
    ut      = y(tt) - x(:,tt)'*b(:,tt);
    Vt      = x(:,tt)'*Pt*x(:,tt) + H;
    Kt      = Tt*Pt*x(:,tt) / Vt;
    Om      = lam*Om + (x(:,tt)*x(:,tt)') + (1-lam)*Om_star;
    Pt      = Om \ In;
    lo(tt)  = (y(tt) - x(:,tt)'*b(:,tt))^2;
    b(:,tt+1) = Tt*b(:,tt) + Kt*ut;
    
end

loss = mean(lo);

end


function  [loss,b,lo] =  rls_EF(para,y,x,pre)

%% initialize and define inputs
[n, T] = size(x);
b  = zeros(n,T);
lo = zeros(T,1);

b(:,1) = 1/n * ones(n,1);
b(1,1) = 0;

lam = para(1);
gam = para(2);

En  = 0.5*eye(n-1) + 0.5*ones(n-1);
Om_ = [1 zeros(1,n-1); zeros(n-1,1), pre*En];
Om  = gam*Om_;
Pt  = Om\eye(n);
H   = lam;
Tt  = eye(n);


%% start recursion
for tt = 1 : T
    
    ut      = y(tt) - x(:,tt)'*b(:,tt);
    Vt      = x(:,tt)'*Pt*x(:,tt) + H;
    Kt      = Tt*Pt*x(:,tt) / Vt;
    Lt      = Tt - Kt*x(:,tt)';
    Pt      = 1/H * (Lt*Pt)';
    lo(tt)  = (y(tt) - x(:,tt)'*b(:,tt))^2;
    b(:,tt+1) = Tt*b(:,tt) + Kt*ut;
    
end

loss = mean(lo);

end


function  [loss,b,lo] = dgas(para,y,x,pre)

%% initialize and define inputs
[n, T] = size(x);
b  = zeros(n,T);
lo = zeros(T,1);

b(:,1) = 1/n*ones(n,1);
b(1,1) = 0;

A1  = speye(n)*para(1)/1e4; A1(1,1) = para(3)/1e4;
B1  = speye(n)*para(2); B1(1,1) = para(4);
om  = (speye(n) - B1)*b(:,1);
kap = para(5);

En = (0.5*eye(n-1) + 0.5*ones(n-1))*pre;
Om = [1 zeros(1,n-1); zeros(n-1,1), En];
In  = eye(n);

Exx     = Om;
xx      = x(:,1)*x(:,1)';
iH      = 0.5*(Exx\In);
sc      = -2 * x(:,1).*(y(1) - x(:,1)'*b(:,1));
b(:,2)  = om + B1*b(:,1) + A1*iH*sc;
lo(1)   = (y(1)-x(:,1)'*b(:,1))^2;


%% start recursion
for tt = 2 : T
    
    Exx     = kap.*Exx+(1-kap).*xx;
    xx      = x(:,tt)*x(:,tt)';
    iH      = 0.5*(Exx\In);
    sc      = -2*x(:,tt).*(y(tt) - x(:,tt)'*b(:,tt));
    lo(tt)  = (y(tt) - x(:,tt)'*b(:,tt))^2;
    b(:,tt+1) = om + B1*b(:,tt) + A1*iH*sc;
    
end

loss = mean(lo);
end


function [L, sig2hat, Sig_star] = ML_est_equ(para, X)

%%
[n, T] = size(X);
Z   = X - mean(X,2);
rho = para(1);

A = sum(diag(Z*Z'));
B = sum(sum(Z).^2);

nu       = rho / (1 + (n-1)*rho);
sig2hat  = (A - nu*B) / (T*n*(1-rho));
Sig_star = sig2hat * ((1-rho)*eye(n) + rho*ones(n,n));

dev = T * (n*log(sig2hat) + (n-1)*log(1-rho) + log(1+(n-1)*rho) + n);
L   = dev/2;

end

%% end of file