-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy patharma_bekk_mvgarchXY4RepeatBayes.m
128 lines (122 loc) · 5.55 KB
/
arma_bekk_mvgarchXY4RepeatBayes.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
function output = arma_bekk_mvgarchXY4RepeatBayes(data, ar_order,bekk_order, indexX, indexY, Nr, Nl, priors, sig, ARMABEKKoptions)
%% QMLE
% exactly, it is mulivariate of ar_arch model at the moment
% model
% X(t) = A(i)X(t-i) + (CC' + B X(t-i) X(t-i)' B')^(1/2) e(t)
% three sequences involved:
% X(t) -- the observations
% E(t) -- the residues of ARMA model
% H(t) -- the conditional variances
% the parameters are estimated by the conditional likelyhood function from
% max(order) to T
% llfn = -0.5ln|H(t)| - 0.5 E(t)'H(t)^(-1)E(t)
% it is the same to minimize -llfn
% input:
% data -- T by k time series
% ar_order -- k by k integer matrix specifying the lags of
% dimension i for dimension j
% bekk_order -- k by k integer matrix specifying the lags of
% dimension i for dimansion j
% Nr, Nl -- number of repeat experiments and the index of the measurements for each
% experiment
% ARMABEKKoptions --
% output:
% output.llfn --- log likelihood function from max(order) to T
% output.parameters -- ARMA, ARCH
% AR(i) = reshape(parameters((i-1)*k^2+i*k^2))
% C = ivech(paramters(k^2*LagAR+1 : k^2*LagAR+1+k(k+1)/2))
% start = k^2*LagAR+1+k(k+1)/2);
% ARCH(i) = reshape(parameters(start+1+(i-1)*k^2:start+1+i*k^2 ))
% output.residues -- E(t)
% output.H -- H(t)
% by now, we are going to use bootstrap to estimate the covariance of the
% parameters
k = size(data,2);
%% estimate the initial parameters by one dimensional model
% fit an one dimensional arma_garch model as the intial paramters for
% amra_bekk
% newA = zeros(k, k, ar_order);
% newC = zeros(k,k);
% newB = zeros(k, k, bekk_order);
% for dim = 1 : k
% spec_fit = garchset('R', ar_order, 'Q', bekk_order, 'Display','off');
% clear tCoeff;
% [tCoeff, tErrors, tLLF] = garchfit(spec_fit, data(:,dim));
% for i = 1 : ar_order
% newA(dim, dim, i) = tCoeff.AR(i)'; % we are using row vector
% end
% newC(dim, dim) = tCoeff.K;
% for i = 1 : bekk_order
% newB(dim,dim, i) = tCoeff.ARCH(i)'; % we are using row vector
% end
% end
% newA = reshape(newA, k*k*ar_order, 1);
[newA, newC, newB] = initialparaforARBEKK(data, k, ar_order, bekk_order);
newC = ivech(newC);
newB = reshape(newB, k,k,bekk_order);
kx = size(indexX,2);
ky = size(indexY,2);
k = kx + ky;
newCx = newC(indexX,indexX);
newCy = newC(indexY,indexY);
newBX = [];
newBY = [];
for i = 1 : bekk_order
newBX(:,:,i) = newB(indexX, :, i);
newBY(:,:,i) = newB(indexY, :, i);
end
startingparameters = [newA;vech(newCx);vech(newCy);reshape(newBX, kx*k*bekk_order,1);reshape(newBY, ky*k*bekk_order, 1)];
Inifeasible = stationary_constraint(startingparameters, ar_order, bekk_order, k, kx, ky);
if Inifeasible > 0
newBX1 = newBX * 0;
newBY1 = newBY * 0;
startingparameters = [newA;vech(newCx);vech(newCy);reshape(newBX1, kx*k*bekk_order,1);reshape(newBY1, ky*k*bekk_order, 1)];
Inifeasible = stationary_constraint(startingparameters, ar_order, bekk_order, k, kx, ky);
if Inifeasible > 0
newA1 = newA * 0;
startingparameters = [newA1;vech(newCx);vech(newCy);reshape(newBX, kx*k*bekk_order,1);reshape(newBY, ky*k*bekk_order, 1)];
Inifeasible = stationary_constraint(startingparameters, ar_order, bekk_order, k, kx, ky);
if Inifeasible > 0
startingparameters = zeros(size(startingparameters));
Inifeasible = stationary_constraint(startingparameters, ar_order, bekk_order, k, kx, ky);
if Inifeasible > 0
'fail to find a feasible initial solution'
end
end
end
end
%% fit the full model by minimizing the llfn
if nargin<=11 || isempty(ARMABEKKoptions)
options.Algorithm = 'interior-point';%'Active-set';
options.LargeScale = 'on';
options.Display='off';
options.Diagnostics='off';
options.TolX=1e-4;
options.TolFun=1e-4;
options.UseParallel = 'Always';
options.MaxFunEvals=5000*length(startingparameters);
options.MaxIter=5000*length(startingparameters);
else
options = ARMABEKKoptions;
end
ObjectFunction = @(parameters)arma_bekk_mvgarch_likelihoodXY4RepeatBayes(parameters, data, ar_order, bekk_order, kx, ky, Nr, Nl, priors, sig);
ConstraintFunction = @(parameters)stationary_constraint(parameters, ar_order, bekk_order, k, kx, ky);
parameters = fmincon(ObjectFunction,startingparameters, [],[],[],[],[],[], ConstraintFunction, options);
% parameters =
% fminunc('arma_bekk_mvgarch_likelihoodXY4Repeat',startingparameters, options, data,ar_order,bekk_order, kx, ky, Nr, Nl);
% ARMABEKKoptions = gaoptimset('PlotFcns',{@gaplotbestf,@gaplotmaxconstr},'Display','iter');
% parameters = ga(FitnessFunction,size(startingparameters,1),[],[],[],[],[],[],ConstraintFunction,ARMABEKKoptions);
%%
[loglikelihood,likelihoods,Ht, errors] = arma_bekk_mvgarch_likelihoodXY4RepeatBayes(parameters,data,ar_order,bekk_order, kx, ky, Nr, Nl, priors, sig);
outputLLF = arma_bekk_mvgarch_likelihoodXY4Repeat4XBayes(parameters, data, ar_order,bekk_order, kx, ky, Nr, Nl, priors, sig, indexX, indexY);
output.LLF_x = -outputLLF.LLF_x;
output.LLF_y = -outputLLF.LLF_y;
output.likelihoods_x = -outputLLF.likelihoods_x;
output.likelihoods_y = -outputLLF.likelihoods_y;
output.Bayes_x = -outputLLF.Bayes_x;
output.Bayes_y = -outputLLF.Bayes_y;
output.Ht = Ht;
output.LLF = -loglikelihood;
output.likelihoods = -likelihoods;
output.parameters = parameters;
output.errors = errors;