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
|
function [bhat,Vq,ACt,uqhat,q2] = clgls(a,yq,Xq,C,qm,mT,qT)
%function [bhat,Vq,ACt,uqhat,q2] = clgls(a,yq,Xq,C,qm,mT,qT)
% bhat = inv[Xq' * inv(Vq) * Xq] * Xq' * inv(Vq) * yq
% Vq = C*A*C'
% ACt = A*C'
% uqhat = yq - Xq * bhat
%
% Written by E.M. Leeper
% Modified by T. Zha, 5/6/97
%
% Copyright (C) 1997-2012 Eric Leeper and Tao Zha
%
% This free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% It is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% If you did not received a copy of the GNU General Public License
% with this software, see <http://www.gnu.org/licenses/>.
%
%*** The following creation of Ahat may be inefficient, T.Z., 5/7/97
Ahat = eye(mT,mT);
i = 1;
arow = zeros(1,mT-1);
while i <= mT-1
arow(1,i) = a^i;
i = i + 1;
end
i = 1;
while i <= mT-1
Ahat(i,i+1:mT) = arow(1,1:mT-i);
i = i + 1;
end
Ahat = Ahat + Ahat' - eye(mT,mT);
%*** GLS to estimate bhat
ACt = Ahat*C';
Vq = C*ACt;
Xqt = Xq';
%CACinv = inv(CAhatC); commented out by T.Z., not necessary
%XqCACinv = Xq' * CACinv; commented out by T.Z., not necessary
bhat = ((Xqt/Vq)*Xq)\(Xqt*(Vq\yq));
uqhat = yq - Xq * bhat;
%*** compute the new quarterly coefficeint "q2"
uqlag = zeros(qT-1,1);
uqlag = uqhat(1:qT-1,1);
uqc = uqhat(2:qT,1);
[u d v]=svd(uqlag,0); %trial
dinv = 1./diag(d); % inv(diag(d))
vd=v.*(ones(size(v,2),1)*diag(d)'); %trial
vdinv=v.*(ones(size(v,2),1)*dinv'); %trial
xtxinv = vdinv*vdinv';
uy = u'*uqc; %trial
xty = vd*uy; %trial
%Bh = xtx\xty; %inv(X'X)*(X'Y), k*m (ncoe*nvar).
q2 = xtxinv*xty; % initial q
|
