/* this will give the likelihood function for each
t to be used with simplex algorithm */
#include <stdio.h>
#include <math.h>
#include "nrutil.h"
#include "factorial.h"
#include "norm_distr.h"
double func_lli(double *p)
{
extern double *index, *rate, *h, *J_stea, *y, *alpha ,
*S , *func_l( double *y_over_h, double *y_over_y ),*temp ;
extern int tot_day, iter ;
int tt, i , k ;
double beta_0, beta_1, beta_2, c_lat , delta, mu_1, gamma_1,
lambda , gamma_2, xi , lli=0.0 ;
double log_k = 0.0 ;
beta_0 = p[1];
beta_1 = p[2];
beta_2 = p[3];
c_lat = p[4];
delta = p[5];
mu_1 = p[6];
gamma_1 = p[7];
lambda = p[8];
gamma_2=sqrt(pow(p[6],2.0)+pow(p[7],2.0));
xi=1+p[8]*pow( sqrt(pow(p[6],2.0)+pow(p[7],2.0)),2.0 );
printf("%s", "gata88");
for( tt=2;tt<=tot_day;tt++ )
{
printf("%s", "gata89");
h[tt]=beta_0 + beta_1*h[tt-1] +
beta_2*h[tt-1]*pow(( J_stea[tt-1]
-lambda*mu_1)/(sqrt(xi))-c_lat,2.0);
if( tt<=tot_day )
{
J_stea[tt]= (log(index[tt]/index[tt-1])-
(rate[tt-1]-h[tt]*(xi)/2.0 ))/sqrt(h[tt]);
y[tt]=sqrt(h[tt])*J_stea[tt];
printf("%s", "gata90");
if ( k>=0 )
{
printf("%s", "gata91");
log_k += log(k);
printf("%s", "gata92");
double temp=l(y[tt]/h[tt],y[tt-1]);
printf("%s", "gata922");
temp +=exp(-lambda+k*log(lambda)-log_k+
log(pgauss(y[tt],k*mu_1*sqrt(h[tt]),h[tt]*(1+k*pow(gamma_1,2)))));
printf("%s", "gata93");
}
}
printf("%s", "gata94");
lli += log(l(y[tt]/h[tt],y[tt-1]));
printf("%s", "gata95");
}
return -(lli);
}
MAIN :
/* this is the main() for the estimation of nlgar
parameters from returns by maximizing the likelihood */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "nrutil.h"
#include "factorial.h"
int N = 8 ;
#define T 1259
int tot_day;
double *index,
*rate,
*h,
/* *J_stea, */
*y,
/* *k_1, */
*l ,
*J_stea,
/* *y_t , */
*S ,
*temp
/* *X_t */
;
FILE *ihandle;
FILE *ohandle,
*ohandle1;
void amoeba(double **,double *,int,double,double(*)(double []),int *);
double func_lli(double *);
double **p,
**e,
*fp,
*x,
lambda ,
gamma_2 ,
*temp_k ,
ftol = 1.0e-8 ;
int nfunk;
int main(void)
{
int i=1,
j,
k=0 ,
counter,
tot_day1
/* N_t */
;
double intr,
mean = 0.0,
sec_mom = 0.0;
double mean_y = 0.0,
mean_h = 0.0,
sec_mom_y = 0.0,
sec_mom_h = 0.0;
double /* mean_X_t=0.0 ,
sec_mom_X =1.0 ,*/
kurt_J_stea = 0.0,
skew_J_stea = 0.0,
sec_mom_skew_J_stea = 0.0,
sec_mom_kurt_J_stea = 0.0;
double var_y,
var_h,
var_kurt_J_stea,
var_skew_J_stea;
double ur,
urn = 0.0,
lvt,
evar=0.0 ;
printf("%s", "gata!");
if(!(ihandle=fopen("c:/test/preturi.txt","rt")))
{
printf("could not open the input file\n");
exit(1);
}
if(!(ohandle=fopen("c:/test/mle.dat","wt")))
{
printf("could not open the first output file\n");
exit(1);
}
if(!(ohandle1=fopen("c:/test/volt.dat","w")))
{
printf("could not open the second output file\n");
exit(1);
}
printf("%s", "gata1");
p=dmatrix(1,N+1,1,N+1);
e=dmatrix(1,N+1,1,N+1);
x=dvector(1,N);
fp=dvector(1,N);
index =dvector(1,T);
rate =dvector(1,T);
while(fscanf(ihandle,"%lf",&index[i])!=EOF)
{
rate[i++] = 0.008807;
}
/* for(j=1;j<N_t;j++)
{
mean_X_t(j)=miu ;
sec_mom_X_T(j)=pow(gamma,2) ;
}
*/
tot_day1=i-1;
for(j=2;j<=tot_day1;j++)
{
sec_mom += pow(log(index[j]/index[j-1]),2.0);
mean += log(index[j]/index[j-1]);
}
tot_day = j-1;
h = dvector(1,tot_day);
J_stea = dvector(1,tot_day);
y = dvector(1,tot_day);
temp = dvector(2,tot_day);
h[1]= sec_mom/(1+lambda*pow(gamma_2,2));
J_stea[1] =0.0;
/*S[1]= 1/sqrt(2*pi) ; */
/* for (j=1 ; j<=Nt, j++)
{
X_t +=
J_stea_t=X_t(0)+sum(1, Nt,X_t(j));
}
for ( N_t=1 , N_t< inf , N_t++ )
{
P(N_t )=exp(-lambda + N_t*log(lambda) - logN_t) ;
}
*/
p[1][1] = 2.0;
p[1][2] = (0.5);
p[1][3] = (400.0);
p[1][4] = (0.2);
p[1][5] = 5.0;
p[1][6] = 2.0 ;
p[1][7] = 1.0 ;
p[1][8] = 3.0 ;
/*
p[1][1] = 1.0;
p[1][2] = (0.6);
p[1][3] = (100.0);
p[1][4] = (0.2);
p[1][5] = 1.0;
*/
printf("%s", "gata2");
for(i=1; i<=N+1; i++)
for(j=1; j<=N+1; j++)
{
e[i][j] = (i==j? 1.0:0.0);
}
for(j=1; j<=N; j++)
x[j] = p[1][j];
fp[1] = func_lli(x);
for(i=2; i<=N+1; i++)
{
for(j=1; j<=N; j++)
{
switch(j)
{
case 1: lambda = 0.5; break;
case 2: lambda = 0.05; break;
case 3: lambda = 5.0; break;
case 4: lambda = 0.5;break;
case 5: lambda = 0.5;
}
x[j] = p[i][j] = p[1][j] + lambda *e[i-1][j];
}
fp[i] = func_lli(x);
}
printf("a\n");
amoeba(p,fp,N,ftol,func_lli,&nfunk);
fprintf(ohandle,"covergence\n");
for(i=1;i<=N+1;i++)
{
for(j=1;j<=N;j++)
{
fprintf(ohandle,"%f ",p[i][j]);
}
fprintf(ohandle,"%f\n",-fp[i]);
}
for(i=2;i<=tot_day;i++)
{
fprintf(ohandle1,"%.2f %.3f\n",index[i],sqrt(h[i+1]*252.0));
printf("%s", "gata3");
/* mean_X_t +=X_t[j]
sec_momX_t +=X_t[j]*X_t[j] ; */
y[i] = sqrt(h[i])*J_stea[i];
/*h[1] = sec_mom/(1+p[1][8]*(pow(p[1][7],2)+pow(p[1][6],2))); */
mean_y +=y[i];
mean_h +=h[i];
sec_mom_y +=y[i]*y[i];
sec_mom_h +=h[i]*h[i];
skew_J_stea += pow(J_stea[i],3.0);
kurt_J_stea += pow(J_stea[i],4.0);
sec_mom_skew_J_stea += pow(J_stea[i],6.0);
sec_mom_kurt_J_stea += pow(J_stea[i],8.0);
}
var_y = (sec_mom_y - tot_day*pow(mean_y/tot_day,2.0))/(tot_day-1);
var_h = (sec_mom_h - tot_day*pow(mean_h/tot_day,2.0))/(tot_day-1);
var_skew_J_stea = (sec_mom_skew_J_stea -
tot_day*pow(skew_J_stea/tot_day,2.0))/(tot_day-1);
var_kurt_J_stea = (sec_mom_kurt_J_stea -
tot_day*pow(kurt_J_stea/tot_day,2.0))/(tot_day-1);
fprintf(ohandle,"variance of y_t = %e\n",var_y);
fprintf(ohandle,"variance of h = %e\n",var_h);
fprintf(ohandle,"skew_J_stea: %lf\n", var_skew_J_stea);
fprintf(ohandle,"kurt_J_stea: %lf\n", var_kurt_J_stea);
printf("%s", "gata4");
ur = p[1][1]+ p[1][2]*(1+pow(p[1][4],2));
fprintf(ohandle1,"unit root condition = %.3f\n",ur);
lvt = p[1][0]/ (1.0 - ur );
fprintf(ohandle,"spot volatility = %.3f\n",sqrt(h[tot_day+1]*252));
fprintf(ohandle,"unconditional mean of volatility = %f\n", lvt );
fclose(ihandle);
fclose(ohandle);
fclose(ohandle1);
printf("%s", "gata5");
system("PAUSE");
}
i have this program ...and it tells at temp that i have a segmentation
fault ...
and i don't know what to do , i looked for the problem ...but ...
