garch.cpp
Go to the documentation of this file.
1 
18 
19 #define WANT_STREAM
20 #define WANT_MATH
21 #define WANT_FSTREAM
22 
23 #include "newmatap.h"
24 #include "newmatio.h"
25 #include "newmatnl.h"
26 
27 #ifdef use_namespace
28 using namespace RBD_LIBRARIES;
29 #endif
30 
31 
32 inline Real square(Real x) { return x*x; }
33 
34 // the class that defines the GARCH log-likelihood
35 
36 class GARCH11_LL : public LL_D_FI
37 {
38  ColumnVector Y; // Y values
39  ColumnVector X; // X values
40  ColumnVector D; // derivatives of loglikelihood
41  SymmetricMatrix D2; // - approximate second derivatives
42  int n; // number of observations
43  Real beta, alpha0, alpha1, beta1;
44  // the parameters
45 
46 public:
47 
48  GARCH11_LL(const ColumnVector& y, const ColumnVector& x)
49  : Y(y), X(x), n(y.nrows()) {}
50  // constructor - load Y and X values
51 
52  void Set(const ColumnVector& p) // set parameter values
53  {
54  para = p;
55  beta = para(1); alpha0 = para(2);
56  alpha1 = para(3); beta1 = para(4);
57  }
58 
59  bool IsValid(); // are parameters valid
60  Real LogLikelihood(); // return the loglikelihood
61  ReturnMatrix Derivatives(); // derivatives of log-likelihood
62  ReturnMatrix FI(); // Fisher Information matrix
63 };
64 
66 { return alpha0>0 && alpha1>0 && beta1>0 && (alpha1+beta1)<1.0; }
67 
69 {
70 // cout << endl << " ";
71 // cout << setw(10) << setprecision(5) << beta;
72 // cout << setw(10) << setprecision(5) << alpha0;
73 // cout << setw(10) << setprecision(5) << alpha1;
74 // cout << setw(10) << setprecision(5) << beta1;
75 // cout << endl;
76  ColumnVector H(n); // residual variances
77  ColumnVector U = Y - X * beta; // the residuals
78  ColumnVector LH(n); // derivative of log-likelihood wrt H
79  // each row corresponds to one observation
80  LH(1)=0;
81  Matrix Hderiv(n,4); // rectangular matrix of derivatives
82  // of H wrt parameters
83  // each row corresponds to one observation
84  // each column to one of the parameters
85 
86  // Regard Y(1) as fixed and don't include in likelihood
87  // then put in an expected value of H(1) in place of actual value
88  // which we don't know. Use
89  // E{H(i)} = alpha0 + alpha1 * E{square(epsilon(i-1))} + beta1 * E{H(i-1)}
90  // and E{square(epsilon(i-1))} = E{H(i-1)} = E{H(i)}
91  Real denom = (1-alpha1-beta1);
92  H(1) = alpha0/denom; // the expected value of H
93  Hderiv(1,1) = 0;
94  Hderiv(1,2) = 1.0 / denom;
95  Hderiv(1,3) = alpha0 / square(denom);
96  Hderiv(1,4) = Hderiv(1,3);
97  Real LL = 0.0; // the log likelihood
98  Real sum1 = 0; // for forming derivative wrt beta
99  Real sum2 = 0; // for forming second derivative wrt beta
100  for (int i=2; i<=n; i++)
101  {
102  Real u1 = U(i-1); Real h1 = H(i-1);
103  Real h = alpha0 + alpha1*square(u1) + beta1*h1; // variance of this obsv.
104  H(i) = h; Real u = U(i);
105  LL += log(h) + square(u) / h; // -2 * log likelihood
106  // Hderiv are derivatives of h with respect to the parameters
107  // need to allow for h1 depending on parameters
108  Hderiv(i,1) = -2*u1*alpha1*X(i-1) + beta1*Hderiv(i-1,1); // beta
109  Hderiv(i,2) = 1 + beta1*Hderiv(i-1,2); // alpha0
110  Hderiv(i,3) = square(u1) + beta1*Hderiv(i-1,3); // alpha1
111  Hderiv(i,4) = h1 + beta1*Hderiv(i-1,4); // beta1
112  LH(i) = -0.5 * (1/h - square(u/h));
113  sum1 += u * X(i)/ h;
114  sum2 += square(X(i)) / h;
115  }
116  D = Hderiv.t()*LH; // derivatives of likelihood wrt parameters
117  D(1) += sum1; // add on deriv wrt beta from square(u) term
118 // cout << setw(10) << setprecision(5) << D << endl;
119 
120  // do minus expected value of second derivatives
121  if (wg) // do only if second derivatives wanted
122  {
123  Hderiv.row(1) = 0.0;
124  Hderiv = H.as_diagonal().i() * Hderiv;
125  D2 << Hderiv.t() * Hderiv; D2 = D2 / 2.0;
126  D2(1,1) += sum2;
127 // cout << setw(10) << setprecision(5) << D2 << endl;
128 // DiagonalMatrix DX; EigenValues(D2,DX);
129 // cout << setw(10) << setprecision(5) << DX << endl;
130 
131  }
132  return -0.5 * LL;
133 }
134 
136 { return D; }
137 
139 {
140  if (!wg) cout << endl << "unexpected call of FI" << endl;
141  return D2;
142 }
143 
144 
145 int my_main()
146 {
147  // get data
148  ifstream fin("garch.dat");
149  if (!fin) { cout << "cannot find garch.dat\n"; exit(1); }
150  int n; fin >> n; // series length
151  // Y contains the dependant variable, X the predictor variable
152  ColumnVector Y(n), X(n);
153  int i;
154  for (i=1; i<=n; i++) fin >> Y(i) >> X(i);
155  cout << "Read " << n << " data points - begin fit\n\n";
156  // now do the fit
157  ColumnVector H(n);
158  GARCH11_LL garch11(Y,X); // loglikehood "object"
159  MLE_D_FI mle_d_fi(garch11,100,0.0001); // mle "object"
160  ColumnVector Para(4); // to hold the parameters
161  Para << 0.0 << 0.1 << 0.1 << 0.1; // starting values
162  // (Should change starting values to a more intelligent formula)
163  mle_d_fi.Fit(Para); // do the fit
164  ColumnVector SE;
165  mle_d_fi.GetStandardErrors(SE);
166  cout << "\n\n";
167  cout << "estimates and standard errors\n";
168  cout << setw(15) << setprecision(5) << (Para | SE) << endl << endl;
169  SymmetricMatrix Corr;
170  mle_d_fi.GetCorrelations(Corr);
171  cout << "correlation matrix\n";
172  cout << setw(10) << setprecision(2) << Corr << endl << endl;
173  cout << "inverse of correlation matrix\n";
174  cout << setw(10) << setprecision(2) << Corr.i() << endl << endl;
175  return 0;
176 }
177 
178 
179 // call my_main() - use this to catch exceptions
180 int main()
181 {
182  Try
183  {
184  return my_main();
185  }
187  {
188  cout << BaseException::what() << "\n";
189  }
190  CatchAll
191  {
192  cout << "\nProgram fails - exception generated\n\n";
193  }
194  return 0;
195 }
196 
197 
#define Try
Definition: myexcept.h:190
#define Catch
Definition: myexcept.h:193
void GetCorrelations(SymmetricMatrix &)
Definition: newmatnl.cpp:256
void Set(const ColumnVector &p)
Definition: garch.cpp:52
ColumnVector Y
Definition: garch.cpp:38
int main()
Definition: garch.cpp:180
double Real
Definition: include.h:307
int my_main()
Definition: garch.cpp:145
#define CatchAll
Definition: myexcept.h:194
int n
Definition: garch.cpp:42
GARCH11_LL(const ColumnVector &y, const ColumnVector &x)
Definition: garch.cpp:48
TransposedMatrix t() const
Definition: newmat6.cpp:320
ReturnMatrix FI()
Definition: garch.cpp:138
DiagedMatrix as_diagonal() const
Definition: newmat6.cpp:339
void Fit(ColumnVector &Parameters)
Definition: newmatnl.cpp:234
bool IsValid()
Definition: garch.cpp:65
The usual rectangular matrix.
Definition: newmat.h:625
InvertedMatrix i() const
Definition: newmat6.cpp:329
ColumnVector X
Definition: garch.cpp:39
static const char * what()
Definition: myexcept.h:93
SymmetricMatrix D2
Definition: garch.cpp:41
Real LogLikelihood()
Definition: garch.cpp:68
ColumnVector D
Definition: garch.cpp:40
void GetStandardErrors(ColumnVector &)
Definition: newmatnl.cpp:253
Column vector.
Definition: newmat.h:1008
GetSubMatrix row(int) const
Definition: submat.cpp:49
Real square(Real x)
Definition: garch.cpp:32
Real beta1
Definition: garch.cpp:43
Symmetric matrix.
Definition: newmat.h:753
ReturnMatrix Derivatives()
Definition: garch.cpp:135


kni
Author(s): Martin Günther
autogenerated on Fri Jan 3 2020 04:01:16