DLResiduals

Purpose

To compute the standardized prediction residuals for a covariance stationary time series.

Old C Code:  hosking.c

This code from the EARMA project can easily be modified to do what we want.

/**********************************************************************/
	input parameters:
	z - vector of length maxlag containing time series observations
	output parameters:
	rllh - likelihood function value
	g - determimant of covariance matrix
	error(maxlag) 
	phiki(maxlag), phikj(maxlag), phikk(maxlag)
***********************************************************************/
#include "stdafx.h"
#include "nrutil.h"
#include <math.h>                                                                                            
double hosking(double *z, double *error, int n, double *c)    	
{	
	double rllh, g, sum;
	int i, j, k;

	VECTOR v;
	VECTOR phi;
	VECTOR phiki;
	VECTOR phikj;

	v = Vector(n);
	phi = Vector(n);
	phiki = Vector(n);
	phikj = Vector(n);
	
	error[0] = z[0];
	v[0] = c[0];
	phi[0] = c[1] / c[0];
	error[1] = z[1] - phi[0] * z[0];
	phiki[0] = phi[0];
	v[1] = v[0] * (1.0 - pow (phi[0], 2));
	g = c[0] * 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 - pow (phi[k - 1], 2));
		g *= v[k];
	}
	sum = 0.0;
	for (j = 0; j< n; j++)
		sum += (pow (error[j], 2)) / v[j]; 
    
	rllh = - 0.5 * ((double)n) * log (sum / ((double)n) ) - 0.5 * log(g); 
	for (i = 0; i < n; i++)
		error[i] /= sqrt(v[i]);

//	free_vector(v);
//	free_vector(phi);     p = 1, q = 0, , phi = 0.6, bfractional = 'n'
//	free_vector(phiki);   p = 1, q = 0, , phi = 0.6, bfractional = 'n'
//	free_vector(phikj);	  p = 1, q = 0, , phi = 0.6, bfractional = 'n'	
	return rllh;

}

Suggested New C Function:  DLResid

I have made the changes I think are needed.  Please polish up.  Remove stdafz.h if possible.  My suggested changes are in magenta.  Please check for memory leaks.

/***************************************************************/
	input parameters:
	z - vector of length n, time series
	n - length of time series
	c - vector of length n containing the autocovariances
	
	output parameters:
	error - vector of length n, prediction residuals
	
	return value:  0, OK
	               1, c is not p.d. 
/****************************************************************/
#include "stdafx.h"
#include "nrutil.h"
#include <math.h>                                                                                            
int DLResid(double *z, double *error, int n, double *c)    	
{	
	double rllh, g, sum;
	int i, j, k;

	VECTOR v;
	VECTOR phi;
	VECTOR phiki;
	VECTOR phikj;

	if (n < 1) return(1)
	v = Vector(n);
	phi = Vector(n);
	phiki = Vector(n);
	phikj = Vector(n);
	
	error[0] = z[0];
	v[0] = c[0];
	if (c[0] <= 0.0) return(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 - pow (phi[0], 2));
	if (v[1] <= 0.0) return(1)
	g = c[0] * 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 - pow (phi[k - 1], 2));
		if (v[k] <= 0.0) return(1)
		}
	for (j = 0; j< n; j++)
		error[j] /= v[j]; 
    
	free_vector(v);
	free_vector(phi);     
	free_vector(phiki);   
	free_vector(phikj);	  
	return(0);

}

Test Example in R

R Function

> DLResiduals
function(r,z)
{
    EPS<-1.110223e-16 # 1+EPS==1, machine epsilon
    n <- length(z)
    if (n!=length(r)) stop("arguments have unequal length")
    error <- numeric(n)
    sigmasq <- numeric(n)
    error[1] <- z[1]
    sigmasq[1] <- r[1]
    phi <- r[2]/r[1]
    error[2] <- z[2] - phi * z[1]
    sigmasqkm1 <- r[1] * (1 - phi^2)
    sigmasq[2] <- sigmasqkm1
    g <- r[1] * sigmasqkm1
    for(k in 2:(n - 1)) {
        if (sigmasqkm1<0 || abs(sigmasqkm1)<EPS) stop("r is not a p.d. sequence")
        phikk <- (r[k + 1] - phi %*% rev(r[2:k]))/sigmasqkm1
        sigmasqk <- sigmasqkm1 * (1 - phikk^2)
        phinew <- phi - phikk * rev(phi)
        phi <- c(phinew, phikk)
        sigmasqkm1 <- sigmasqk
        g <- g * sigmasqk
        error[k + 1] <- z[k + 1] - crossprod(phi, rev(z[1:k]))
        sigmasq[k + 1] <- sigmasqk
    }
    error/sqrt(sigmasq)
}

Example

> z<-1:5
> DLResiduals(1/z, z)
[1] 1.000000 1.732051 2.323790 2.861920 3.368777