varve = ts(scan("http://www.stat.ualberta.ca/~wiens/stat479/S&Sdatasets/varve.dat"))
varve = ts(scan("varve.dat"))

## Set the CRAN mirror to Canada (BC) (Menu:  Packages -> Set CRAN mirror)
## Install package fracdiff (Packages -> Install package)
## Load it: Packages -> Load package
library(fracdiff)

lvarve = log(varve) - mean(log(varve)) # Centre the logs
par(mfrow=c(3,1))

acf(lvarve, 50)

## Estimate d:
varve.fd = fracdiff(lvarve, nar=0, nma=0)
summary(varve.fd)

names(varve.fd)

d <- varve.fd$d	
cat("estimated d =", d, "\n")
cat("std.err of d.hat = ", varve.fd$stderror.dpq, "\n")

## Compute M residuals: 
M = 500
## First compute the pi's (P_j = pi_(j-1)):
P = vector(length=M)
P[1] = 1
for(j in 1:(M-1)) P[j+1] = (j-1-d)*P[j]/j
## Now the residuals:
resid = vector(length = M)
for(t in 1:M) resid[t] = P[1:t]%*%lvarve[t:1]
resid = resid[101:M]  #Remove the first 100


plot(resid)
acf.varve = acf(resid, lag.max = 50)