%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% QMLE estimation of the excess (with respect to the risk-free rate) return process:
% The conditional mean is driven by the dividend yield (to evaluate predictability)
% and the error terms are given by a GARCH(1,1) process
Copyright By PowCoder代写 加微信 powcoder
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Clean the environment
close all;
clear all;
% Upload the data
data = readtable(‘predictability.xls’);
ret_m = log(1+data.retdny); % Continuously-compounded market returns
ret_f = log(1+data.tbill); % Continuously-compounded riskless returns
dp = data.dyny; % The divinded yield
r = ret_m – ret_f; % Excess continuously-compounded returns (risk premia)
T = length(r);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Let us first estimate the conditional mean by least squares:
% This is useful for two reasons:
% (1) It is helpful to have a benchmark.
% (2) You need starting values for the QMLE procedure. It makes therefore
% sense to use the OLS values as starting values.
% We will also compute OLS standard errors: again, to have a benchmark.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
X = [ones(T-1,1),dp(1:end-1)];
beta = inv(X’*X)*(X’*r(2:end));
OLSres = r(2:end) – X*beta;
sigma2 = (OLSres’*OLSres)/(length(OLSres)-2);
Var_beta = sigma2*inv(X’*X);
OLStstat = beta./sqrt(diag(Var_beta));
OLStstat = [OLStstat;0;0;0];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Always plot the data (when you can):
% Do we need a time-varying conditional mean?
% Do we need a time-varying conditional variance?
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
plot(r(2:end));
plot(beta(1)+beta(2)*dp(1:end-1));
title(‘Excess continuously-compounded returns and their conditional mean’)
ylabel(‘returns and conditional expected returns’)
xlabel(‘time’)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Let us now maximize the log-likelihood and find the parameter estimates
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
parameter_names = [“alpha” “beta” “mu” “delta” “phi”]’;
% Initial conditions (I am just choosing reasonable values, but you can of course change them and see what happens)
theta_initial=[beta(1),beta(2),0,0.5,0.5];
options = optimset(‘Display’,’none’,’TolFun’,1e-20,’TolX’,1e-20);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Let us now do inference: we need the asymptotic variance/covariance
% We will compute it in its “sandwich” format.
% Recall: QMLE requires both the gradient and the hessian
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sum4 = zeros(5,5);
sum5 = zeros(5,5);
h = var(r)*ones(T-1,1); % We are initializing the conditional variance at the unconditional level
grad1_h = zeros(T-1,1); % We are initializing the first derivatives at zero
grad2_h = zeros(T-1,1);
grad3_h = zeros(T-1,1);
hes12_h = zeros(T-1,1); % We are initializing the second derivatives at zero
hes22_h = zeros(T-1,1);
hes23_h = zeros(T-1,1);
h(t) = theta(3)+theta(4)*h(t-1)+theta(5)*(r(t-1)-(theta(1)+theta(2)*dp(t-2)))^2;
% Compute the gradient: this is the first piece (2 rows) having to do with the conditional mean
grad1 = [(r(t) – (theta(1)+theta(2)*dp(t-1)))/h(t); (r(t) – (theta(1)+theta(2)*dp(t-1)))*dp(t-1)/h(t)]; % First piece
% Compute the gradient: this is the second piece (3 rows) having to do with the conditional variance
% First, we have to have a recursion over the first derivatives of h with respect to the parameters
grad1_h(t) = 1 + theta(4)*grad1_h(t-1);
grad2_h(t) = h(t-1) + theta(4)*grad2_h(t-1);
grad3_h(t) = (r(t-1)-(theta(1)+theta(2)*dp(t-2)))*(r(t-1)-(theta(1)+theta(2)*dp(t-2))) + theta(4)*grad3_h(t-1);
grad_h = [grad1_h(t); grad2_h(t); grad3_h(t)];
grad = -(1/2)*1/h(t)+(1/2)*(r(t)-(theta(1)+theta(2)*dp(t-1)))*(r(t)-(theta(1)+theta(2)*dp(t-1)))/(h(t)^2);
grad2 = grad.*grad_h; % Second piece
gradient = [grad1;grad2]; % Full gradient – first piece on top of the second piece
sum4 = sum4+gradient*gradient’; % Loop to compute Omega_0
% Compute the hessian:
% First, we need to have a recursion over some of the second derivatives of h with respect to the parameters
v = (1/2)*(1/(h(t)^2)) – ((r(t)-(theta(1)+theta(2)*dp(t-1)))^2)/(h(t)^3);
hes12_h(t) = grad1_h(t-1) + theta(4)*hes12_h(t-1);
hes22_h(t) = 2*grad2_h(t-1) + theta(4)*hes22_h(t-1);
hes23_h(t) = grad3_h(t-1) + theta(4)*hes23_h(t-1);
hessian = [-1/h(t), -dp(t-1)/h(t), 0, 0, 0; …
-dp(t-1)/h(t), -dp(t-1)^2/h(t),0,0,0; …
0, 0, v*(grad1_h(t))^2, v*(grad1_h(t))*(grad2_h(t)) + grad*hes12_h(t), v*(grad1_h(t))*(grad3_h(t)); …
0, 0, v*(grad2_h(t))*(grad1_h(t))+grad*hes12_h(t), v*(grad2_h(t))^2 + grad*hes22_h(t), v*(grad2_h(t))*(grad3_h(t))+grad*hes23_h(t); …
0, 0, v*(grad3_h(t))*(grad1_h(t)), v*(grad3_h(t))*(grad2_h(t))+grad*hes23_h(t),v*(grad3_h(t))^2];
sum5=sum5+hessian;
asy_varMLE = (1/T)*inv((1/T)*sum4); % Asymptotic variance under correct specification
t_statsMLE = theta’./sqrt(diag(asy_varMLE));
asy_varQMLE = (1/T)*inv(-(1/T)*sum5)*((1/T)*sum4)*inv(-(1/T)*sum5); % QMLE asymptotic variance
t_statsQMLE = theta’./sqrt(diag(asy_varQMLE));
outcome = table(parameter_names,theta’,t_statsMLE, t_statsQMLE, OLStstat, …
‘VariableNames’,{‘parameters’ ‘Estimates_QMLE’ ‘t_statistics_MLE’ ‘t_statistics_QML’ ‘t_statistics_OLS’});
disp(outcome);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This is the function which delivers the log-likelihood given an assumption
% of normality on the error terms
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y=garchloglik(theta,r,dp)
h1 = var(r); % Initial value of the conditional variance: we use the unconditional
% empirical variance from the data
T = length(r);
alpha = theta(1);
beta = theta(2);
mu = theta(3);
delta = theta(4);
phi = theta(5);
sum1=-T/2*log(2*pi);
sum2=-0.5*log(h1);
sum3= -(r(2)-(alpha+beta*dp(1)))*(r(2)-(alpha+beta*dp(1)))/(2*h1);
h = mu+delta*h+phi*(r(t-1)-(alpha+beta*dp(t-2)))*(r(t-1)-(alpha+beta*dp(t-2)));
sum2=sum2-0.5*log(h);
stand_res=(r(t)-(alpha+beta*dp(t-1)))*(r(t)-(alpha+beta*dp(t-1)))/(2*h);
sum3=sum3-stand_res;
y= -(1/T)*(sum1+sum2+sum3); % Note that I changed the sign of the likelihood. We are maximizing, not minimizing.
% The function garchloglik is however fed into a
% minsearch function
程序代写 CS代考 加微信: powcoder QQ: 1823890830 Email: powcoder@163.com