https://github.com/cran/ltsa
Tip revision: 7e0165a782fedb4fc5aaa42677ea20ad995930a3 authored by A.I. McLeod on 30 November 2010, 19:46:10 UTC
version 1.4
version 1.4
Tip revision: 7e0165a
DLResid.c
/**************************************************************************************
DLResid
input parameters:
z - vector of length n, time series
nR - length of time series
c - vector of length n containing the autocovariances
EPSL - machine epsilon in double precision
stanQ - =0, raw residuals otherwise standardized
output parameters:
error - vector of length n, prediction residuals
LogLike - loglikelihood
ifault: 0, OK; 1, c is not p.d.
***************************************************************************************/
#include "nrutil.h"
#include <math.h>
void DLResid(double *z,double *error, int *nR,double *c, double *LogLike, double *EPSL, int *stanQ, int *fault)
{
double logg,sum;
int i,j,k,n;
double EPS;
VECTOR v,phi,phiki,phikj;
*fault = 0;
n = *nR;
EPS = *EPSL;
if (n < 1) *fault = 1;
v = Vector(n);
phi = Vector(n);
phiki = Vector(n);
phikj = Vector(n);
error[0] = z[0];
v[0] = c[0];
if (c[0] <= EPS) *fault = 1;
phi[0] = c[1] / c[0];
error[1] = z[1] - phi[0] * z[0];
phiki[0] = phi[0];
v[1] = v[0] * (1.0 - phi[0]*phi[0]);
if (v[1] <= EPS) *fault = 1;
logg = log(c[0]) + log(v[1]);
for (k = 2; k < n; k++)
{
sum = 0.0;
for (i = 1; i < k; i++)
sum += phiki[i - 1] * c[k - i];
phi[k - 1] = (c[k] - sum) / v[k - 1];
for (j = 1; j< k; j++)
{
phikj[j - 1] = phiki[j - 1] - phi[k - 1] * phiki[k - j - 1];
}
phikj[k - 1] = phi[k - 1];
sum = 0.0;
for (j = 1; j<=k; j++)
{
sum += phikj[j - 1] * z[k - j];
phiki[j - 1] = phikj[j - 1];
}
error[k] = z[k] - sum;
v[k] = v[k - 1] * (1.0 - phi[k - 1]*phi[k - 1]);
logg += log(v[k]);
if (v[k] <= EPS) *fault = 1;
}
for (j = 0; j< n; j++)
error[j] /= sqrt(v[j]);
sum = 0.0;
for (j = 0; j< n; j++)
sum += error[j]*error[j];
*LogLike = - 0.5 * ((double)n) * log (sum / ((double)n) ) - 0.5 * logg;
if ( !*stanQ ) {
for (j = 0; j< n; j++)
error[j] *= sqrt(v[j]);
}
free_vector(v);
free_vector(phi);
free_vector(phiki);
free_vector(phikj);
return;
}