generate_fogm_cov.c

最大似然估计算法

#include "timeseries.h"

void generate_FOGM_cov(double *time, int n_pts, double Delta_t, double FOGM, double beta, double *C)
{

	int j, k;

	double sec_per_year, fs, alpha, gamma, a, b, t_prev;
	double sigma_W;

	double *dt, *temp1, *work1;

	alpha = 1.0;
	gamma = 0.0;
	sec_per_year = 365.24219*24.0*3600.0;

	beta /= sec_per_year;

	dt    = (double *)calloc((size_t)(n_pts-1),     sizeof(double));
	temp1 = (double *)calloc((size_t)(n_pts),       sizeof(double));
	work1 = (double *)calloc((size_t)(n_pts*n_pts), sizeof(double));

	for (j = 0; j < n_pts-1; j++) dt[j] = time[j+1] - time[j];

	fs = 1.0 / Delta_t / sec_per_year;

	temp1[0] = sqrt(Delta_t);
	for (j = 1; j < n_pts; j++) {
		temp1[j] = sqrt(dt[j-1]);
	}

	a = -beta * sec_per_year;
	for (j = 0; j < n_pts; j++) {
		for (k = 0; k < n_pts; k++) {
			if (j >= k) work1[j + k * n_pts] = exp(a * (time[j] - time[k]) );
			else work1[j + k * n_pts] = 0.0;
		}
	}


	sigma_W = FOGM / beta * sqrt(fs/sec_per_year);

	a = exp(-beta/fs) * FOGM;
	b = exp(-fs/beta) * sigma_W;

	for (j = 0; j < n_pts; j++) temp1[j] = a * temp1[j] + b;

	for (j = 0; j < n_pts; j++) {
		for (k = 0; k < n_pts; k++) {
			work1[j + k * n_pts] *= temp1[k];
		}
	}

	dgemm_("N","T",&n_pts,&n_pts,&n_pts,&alpha,work1,&n_pts,work1,&n_pts,&gamma,C,&n_pts);

	free(dt);
	free(temp1);
	free(work1);
	
}