Go to the documentation of this file.00001
00002
00003
00018
00019 #define WANT_STREAM
00020 #define WANT_MATH
00021 #define WANT_FSTREAM
00022
00023 #include "newmatap.h"
00024 #include "newmatio.h"
00025 #include "newmatnl.h"
00026
00027 #ifdef use_namespace
00028 using namespace RBD_LIBRARIES;
00029 #endif
00030
00031
00032 inline Real square(Real x) { return x*x; }
00033
00034
00035
00036 class GARCH11_LL : public LL_D_FI
00037 {
00038 ColumnVector Y;
00039 ColumnVector X;
00040 ColumnVector D;
00041 SymmetricMatrix D2;
00042 int n;
00043 Real beta, alpha0, alpha1, beta1;
00044
00045
00046 public:
00047
00048 GARCH11_LL(const ColumnVector& y, const ColumnVector& x)
00049 : Y(y), X(x), n(y.nrows()) {}
00050
00051
00052 void Set(const ColumnVector& p)
00053 {
00054 para = p;
00055 beta = para(1); alpha0 = para(2);
00056 alpha1 = para(3); beta1 = para(4);
00057 }
00058
00059 bool IsValid();
00060 Real LogLikelihood();
00061 ReturnMatrix Derivatives();
00062 ReturnMatrix FI();
00063 };
00064
00065 bool GARCH11_LL::IsValid()
00066 { return alpha0>0 && alpha1>0 && beta1>0 && (alpha1+beta1)<1.0; }
00067
00068 Real GARCH11_LL::LogLikelihood()
00069 {
00070
00071
00072
00073
00074
00075
00076 ColumnVector H(n);
00077 ColumnVector U = Y - X * beta;
00078 ColumnVector LH(n);
00079
00080 LH(1)=0;
00081 Matrix Hderiv(n,4);
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091 Real denom = (1-alpha1-beta1);
00092 H(1) = alpha0/denom;
00093 Hderiv(1,1) = 0;
00094 Hderiv(1,2) = 1.0 / denom;
00095 Hderiv(1,3) = alpha0 / square(denom);
00096 Hderiv(1,4) = Hderiv(1,3);
00097 Real LL = 0.0;
00098 Real sum1 = 0;
00099 Real sum2 = 0;
00100 for (int i=2; i<=n; i++)
00101 {
00102 Real u1 = U(i-1); Real h1 = H(i-1);
00103 Real h = alpha0 + alpha1*square(u1) + beta1*h1;
00104 H(i) = h; Real u = U(i);
00105 LL += log(h) + square(u) / h;
00106
00107
00108 Hderiv(i,1) = -2*u1*alpha1*X(i-1) + beta1*Hderiv(i-1,1);
00109 Hderiv(i,2) = 1 + beta1*Hderiv(i-1,2);
00110 Hderiv(i,3) = square(u1) + beta1*Hderiv(i-1,3);
00111 Hderiv(i,4) = h1 + beta1*Hderiv(i-1,4);
00112 LH(i) = -0.5 * (1/h - square(u/h));
00113 sum1 += u * X(i)/ h;
00114 sum2 += square(X(i)) / h;
00115 }
00116 D = Hderiv.t()*LH;
00117 D(1) += sum1;
00118
00119
00120
00121 if (wg)
00122 {
00123 Hderiv.row(1) = 0.0;
00124 Hderiv = H.as_diagonal().i() * Hderiv;
00125 D2 << Hderiv.t() * Hderiv; D2 = D2 / 2.0;
00126 D2(1,1) += sum2;
00127
00128
00129
00130
00131 }
00132 return -0.5 * LL;
00133 }
00134
00135 ReturnMatrix GARCH11_LL::Derivatives()
00136 { return D; }
00137
00138 ReturnMatrix GARCH11_LL::FI()
00139 {
00140 if (!wg) cout << endl << "unexpected call of FI" << endl;
00141 return D2;
00142 }
00143
00144
00145 int my_main()
00146 {
00147
00148 ifstream fin("garch.dat");
00149 if (!fin) { cout << "cannot find garch.dat\n"; exit(1); }
00150 int n; fin >> n;
00151
00152 ColumnVector Y(n), X(n);
00153 int i;
00154 for (i=1; i<=n; i++) fin >> Y(i) >> X(i);
00155 cout << "Read " << n << " data points - begin fit\n\n";
00156
00157 ColumnVector H(n);
00158 GARCH11_LL garch11(Y,X);
00159 MLE_D_FI mle_d_fi(garch11,100,0.0001);
00160 ColumnVector Para(4);
00161 Para << 0.0 << 0.1 << 0.1 << 0.1;
00162
00163 mle_d_fi.Fit(Para);
00164 ColumnVector SE;
00165 mle_d_fi.GetStandardErrors(SE);
00166 cout << "\n\n";
00167 cout << "estimates and standard errors\n";
00168 cout << setw(15) << setprecision(5) << (Para | SE) << endl << endl;
00169 SymmetricMatrix Corr;
00170 mle_d_fi.GetCorrelations(Corr);
00171 cout << "correlation matrix\n";
00172 cout << setw(10) << setprecision(2) << Corr << endl << endl;
00173 cout << "inverse of correlation matrix\n";
00174 cout << setw(10) << setprecision(2) << Corr.i() << endl << endl;
00175 return 0;
00176 }
00177
00178
00179
00180 int main()
00181 {
00182 Try
00183 {
00184 return my_main();
00185 }
00186 Catch(BaseException)
00187 {
00188 cout << BaseException::what() << "\n";
00189 }
00190 CatchAll
00191 {
00192 cout << "\nProgram fails - exception generated\n\n";
00193 }
00194 return 0;
00195 }
00196
00197