#This is the S-Plus function.
#Written by Vadim teverovsky

Sdurbin.FGN <- function(num,h,sigma)

# FUNCTION which produces a FGN sequence by calling durbin_FGN.c.
# It produces a normal sequence, for the innovations, and then lets C do
# most of the work.  Vadim Teverovsky 1994.
#	Inputs:
#		num	-- length of desired series.
#		h	-- Hurst parameter.
#		sigma	-- the standard deviation of the innovations used.
#	Output:
#		out 	-- FGN vector of length num.
#

{
	if(!is.loaded(C.symbol("durbin_FGN.o")))
	dyn.load("durbin_FGN.o")
	normal <- rnorm(num+1,0,1)
	out<-rep(0,num+1)
	.C("durbinFGN",
		num = as.integer(num),
		h = as.double(h),
		sigma = as.double(sigma),
		normal = as.double(normal),
		out = as.double(out))$out[2:(num+1)]
}



#####This is the beginning of the C function.  The object file needs to be
#produced for the S routine above.

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

/*Function to simulate a FGN  sequence, using the Durbin-Levinson 
coefficients, and working in conjunction with Sdurbin.FGN.S, an S-Plus 
routine.  */

void durbinFGN(n,h,sigma,normal,output)

long *n;
double *h, *sigma, *normal, *output;

/* n  -- length of desired sequence.
   d  -- coefficient
   sigma -- standard deviation of the innovations. 
   normal  -- the sequence of innovations.
   output -- the sequence being output.*/

{
  long i,j;
  double sigma2, *acov, *vee, *phi1, *phi2, clock1, clock2;
  acov = (double *) S_alloc(*n+1,sizeof(double));
  vee  = (double *) S_alloc(*n+1,sizeof(double));
  phi1  = (double *) S_alloc(*n+1,sizeof(double));
  phi2  = (double *) S_alloc(*n+1,sizeof(double));
  
/*    All the vectors necessary are allocated dynamically.*/

  sigma2 = *sigma * (*sigma)/2;

/*determining the autocovariance function for this particular h.*/

  for (i = 0 ; i <= *n; i++)
    acov[i] = sigma2*(pow((double)(i+1),2*(*h)) - 2*pow((double)i,2*(*h)) + 
		      pow((double)abs(i-1),2*(*h)));

/* Determining the Durbin-Levinson coefficients and the output vector, recursively*/

  phi1[1] = acov[1]/acov[0];
  phi2[1] = phi1[1];
  vee[0] = acov[0];
  vee[1] = vee[0]*(1 - phi1[1]*phi1[1]);
  output[1] = sqrt(vee[0])*normal[1];

  for (i=2; i <= *n; i++){
    phi1[i] = acov[i];

    for (j=1; j <= i-1; j++)
      phi1[i] -= phi2[j]*acov[i-j];
    phi1[i] = phi1[i]/vee[i-1];

    vee[i] = vee[i-1]*(1-phi1[i]*phi1[i]);
    output[i] = sqrt(vee[i-1])*normal[i];
    for (j=1; j <= i-1; j++){
      phi1[j] = phi2[j]-phi1[i]*phi2[i-j];

      output[i] += phi2[j]*output[i-j];
    }
    
    for (j=1; j <= i; j++)
      phi2[j] =  phi1[j];
  }

  free(acov);
  free(vee);
  free(phi1);
  free(phi2);
}