Fractional Differentiation (R) — the fracdiff package cross-check¶
diffseries for the transform, fdGPH/fracdiff to estimate the memory parameter d¶
The R companion to the from-scratch fractional differentiation. Where Python has no maintained package for this (and mlfinlab is commercial), R's CRAN fracdiff package (Fraley, Leisch, Maechler; Hosking's method) is the authoritative reference: diffseries(x, d) applies the fractional-difference operator, and fdGPH / fdSperio / fracdiff estimate the long-memory parameter $d$ from the data. We reproduce the stationarity-vs-memory frontier with the package — confirming the Python from-scratch result — and add what the package uniquely offers: a data-driven estimate of $d$. Same S&P 500 log-price. This is the package half of the "from-scratch + package" pairing.
options(repr.plot.width=9, repr.plot.height=4.6, warn=-1)
.libPaths(c("C:/Users/user/R/win-library/4.6", .libPaths()))
for(p in c("tseries")) if(!requireNamespace(p,quietly=TRUE)) install.packages(p, repos="https://cloud.r-project.org", lib="C:/Users/user/R/win-library/4.6", quiet=TRUE)
suppressMessages({library(fracdiff); library(tseries)})
r <- read.csv("spx_rv_ret.csv")$ret/100
lp <- cumsum(r) # S&P log price
cat("fracdiff", as.character(packageVersion("fracdiff")), "| log-price length", length(lp), "\n")
Registered S3 method overwritten by 'quantmod': method from as.zoo.data.frame zoo
fracdiff 1.5.4 | log-price length 3459
1. diffseries — the frontier, from the package¶
fracdiff::diffseries(x, d) fractionally differences the series by order $d$. We sweep $d\in[0,1]$, and at each value run the augmented Dickey-Fuller test (tseries::adf.test) and the correlation with the original price — the same frontier the Python notebook built by hand.
The frontier has the same shape, and the minimum-$d$ answer does not match: the package route selects a noticeably smaller $d$ than the from-scratch one. That is worth understanding rather than papering over, because neither implementation is wrong. diffseries applies the operator over the whole available history (an expanding window), while the Python version truncates the weights at a threshold and uses a fixed-width window — López de Prado's own recommendation, since it keeps the transform's memory length constant over the sample. The two also call different ADF tests: tseries::adf.test fixes the lag order at $\lfloor(n-1)^{1/3}\rfloor$, while statsmodels selects it by AIC.
So two defensible implementations of the same recipe disagree about the answer by several grid steps. Section 2 shows that the disagreement does not actually matter, because both are below the order at which the series becomes stationary at all.
ds <- seq(0,1,by=0.05); adfp <- numeric(length(ds)); corr <- numeric(length(ds))
for(i in seq_along(ds)){
fd <- diffseries(lp, ds[i])
adfp[i] <- suppressWarnings(adf.test(fd)$p.value)
corr[i] <- cor(fd, lp)
}
mind <- ds[which(adfp<0.05)[1]]
par(mar=c(4,4,3,4))
plot(ds, adfp, type="b", pch=19, col="#c53030", lwd=2, xlab="fractional order d", ylab="ADF p-value", main=sprintf("fracdiff::diffseries frontier -- min stationary d ~ %.2f", mind))
abline(h=0.05, lty=3, col="#c53030"); abline(v=mind, lty=2, col="#2f855a")
par(new=TRUE); plot(ds, corr, type="b", pch=15, col="#2b6cb0", lwd=2, axes=FALSE, xlab="", ylab=""); axis(4); mtext("correlation with price (memory)", 4, 3, col="#2b6cb0")
cat(sprintf("log-price ADF p %.2f (non-stationary); min-d %.2f: ADF p<0.05 AND memory %.2f; returns d=1 memory %.2f\n",
adfp[1], mind, corr[which(adfp<0.05)[1]], corr[length(corr)]))
cat(sprintf("\nThe CRAN package reproduces the SHAPE of the frontier, and picks d = %.2f where the from-scratch\n", mind))
cat("fixed-width version picked 0.35. Same recipe, different implementations of it: diffseries uses an expanding\n")
cat("window over all history, the Python version a truncated fixed-width one, and the two ADF routines choose their\n")
cat("lag orders differently. A rule whose answer moves this much between two correct implementations is a rule\n")
cat("carrying more weight than it can bear -- which is the subject of the next section.\n")
log-price ADF p 0.63 (non-stationary); min-d 0.25: ADF p<0.05 AND memory 0.91; returns d=1 memory 0.05
The CRAN package reproduces the SHAPE of the frontier, and picks d = 0.25 where the from-scratch
fixed-width version picked 0.35. Same recipe, different implementations of it: diffseries uses an expanding
window over all history, the Python version a truncated fixed-width one, and the two ADF routines choose their
lag orders differently. A rule whose answer moves this much between two correct implementations is a rule
carrying more weight than it can bear -- which is the subject of the next section.
2. Estimating d from the data — the integration order vs the minimum-d¶
This is what the package offers that a hand-rolled operator does not: it can estimate the memory parameter rather than assume it. fdGPH runs the Geweke-Porter-Hudak log-periodogram regression, fdSperio a smoothed variant, and fracdiff() fits an ARFIMA(0,$d$,0) by maximum likelihood.
Applied to the log price these agree it is a near-random-walk. Applied to the differenced series they answer the question the ADF test could not, and this is the check that settles the section. A fractionally integrated $I(\delta)$ process is stationary exactly when $\delta<0.5$ — no test required, it is a property of the process. So estimating $\delta$ of each candidate series says directly whether the minimum-$d$ rule chose an order that works.
One caution on reading the three estimators. fracdiff()'s maximum likelihood is constrained to $d<0.5$ because ARFIMA is only defined as a stationary model there, so on an $I(1)$ series it returns the boundary rather than an estimate. A value of exactly 0.500 is the optimiser pinned against its constraint, not a measurement, and it should not be averaged in with the other two as though it were a third opinion.
gph <- fdGPH(lp)$d
sper <- fdSperio(lp)$d
mle <- fracdiff(lp, nar=0, nma=0)$d
cat(sprintf("Integration order of the LOG PRICE: fdGPH %.3f | fdSperio %.3f | ARFIMA-MLE %.3f\n", gph, sper, mle))
cat(sprintf(" GPH and Sperio agree the price is ~I(1), a near-random-walk. The MLE returns exactly %.3f because that\n", mle))
cat(" is its constraint boundary -- ARFIMA is only defined for d<0.5, so on an I(1) series it reports the edge of\n")
cat(" the parameter space rather than an estimate. Two opinions here, not three.\n\n")
cat("Now the question the ADF test could not answer -- what is the order of the DIFFERENCED series?\n")
cat("An I(delta) process is stationary exactly when delta < 0.5.\n\n")
cat(sprintf(" %-22s %12s %14s %13s\n", "series", "fdGPH order", "stationary?", "memory kept"))
cat(sprintf(" %-22s %12.3f %14s %13.2f\n", "log price (d=0)", gph, "NO", 1.00))
for(dd in c(0.25,0.30,0.35,0.40,0.45,0.50,0.55,0.60)){
fdx <- diffseries(lp, dd); g <- fdGPH(fdx)$d
cat(sprintf(" %-22s %12.3f %14s %13.2f\n", sprintf("diffseries d=%.2f",dd), g,
ifelse(g<0.5,"yes","NO"), cor(fdx, lp)))
}
g1 <- fdGPH(diff(lp))$d
cat(sprintf(" %-22s %12.3f %14s %13.2f\n", "returns (d=1)", g1, ifelse(g1<0.5,"yes","NO"), cor(c(NA,diff(lp))[-1], lp[-1])))
cat(sprintf("\n The order falls roughly one-for-one with d, as the theory requires. The series crosses below 0.5 at\n"))
cat(sprintf(" around d = 0.50 -- well ABOVE the d = %.2f the ADF frontier selected here, and above the 0.35 the\n", mind))
cat(" Python notebook selected. Both minimum-d answers are too small: those series pass an ADF test and are\n")
cat(" still non-stationary. This is the package earning its place -- fdGPH answers directly what adf.test\n")
cat(" could only guess at, and it agrees with the from-scratch GPH in the Python companion to three decimals.\n")
cat("\n The method itself is untouched. At d=0.50 the series is genuinely stationary and still retains far more\n")
cat(" of the price level than returns do. Only the selection rule was wrong.\n")
# show the chosen fracdiff series vs price vs returns, at the order that is actually stationary
dhat <- 0.50
fd <- diffseries(lp, dhat)
par(mfrow=c(3,1), mar=c(2,4,2,1))
plot(lp, type="l", col="#a0aec0", ylab="log price", main="d=0: log price (non-stationary)")
plot(fd, type="l", col="#2f855a", ylab=sprintf("fracdiff d=%.2f",dhat), main=sprintf("d=%.2f: stationary by integration order, memory retained (diffseries)",dhat)); abline(h=mean(fd),col="black",lwd=.5)
plot(c(NA,diff(lp)), type="l", col="#c53030", ylab="returns", main="d=1: returns (memory destroyed)"); abline(h=0,col="black",lwd=.5)
par(mfrow=c(1,1))
Integration order of the LOG PRICE: fdGPH 0.926 | fdSperio 0.992 | ARFIMA-MLE 0.500
GPH and Sperio agree the price is ~I(1), a near-random-walk. The MLE returns exactly 0.500 because that
is its constraint boundary -- ARFIMA is only defined for d<0.5, so on an I(1) series it reports the edge of
the parameter space rather than an estimate. Two opinions here, not three.
Now the question the ADF test could not answer -- what is the order of the DIFFERENCED series?
An I(delta) process is stationary exactly when delta < 0.5.
series fdGPH order stationary? memory kept
log price (d=0) 0.926 NO 1.00
diffseries d=0.25 0.691 NO 0.91 diffseries d=0.30 0.657 NO 0.86 diffseries d=0.35 0.619 NO 0.80 diffseries d=0.40 0.578 NO 0.73 diffseries d=0.45 0.541 NO 0.65 diffseries d=0.50 0.498 yes 0.56 diffseries d=0.55 0.452 yes 0.47 diffseries d=0.60 0.406 yes 0.38
returns (d=1) -0.059 yes 0.05
The order falls roughly one-for-one with d, as the theory requires. The series crosses below 0.5 at
around d = 0.50 -- well ABOVE the d = 0.25 the ADF frontier selected here, and above the 0.35 the
Python notebook selected. Both minimum-d answers are too small: those series pass an ADF test and are
still non-stationary. This is the package earning its place -- fdGPH answers directly what adf.test
could only guess at, and it agrees with the from-scratch GPH in the Python companion to three decimals.
The method itself is untouched. At d=0.50 the series is genuinely stationary and still retains far more
of the price level than returns do. Only the selection rule was wrong.
3. Summary¶
The CRAN fracdiff package reproduces the shape of the frontier and then does something the from-scratch operator cannot: it measures the integration order instead of inferring it from a hypothesis test.
That measurement is what corrects the chapter. diffseries and the hand-built fixed-width operator disagree about the minimum ADF-stationary $d$ — the package route picks a smaller one — and both answers turn out to be too small anyway. Running fdGPH on the differenced series shows the order crossing below the $0.5$ stationarity boundary only at around $d=0.50$; at the ADF-selected orders the result still measures well above it. Those series pass a unit-root test and are not stationary.
The package's own estimators need reading with care too: fdGPH and fdSperio put the log price near $I(1)$, while fracdiff()'s maximum likelihood returns exactly 0.500 because ARFIMA is only defined for $d<0.5$ and the optimiser is pinned to its boundary. That is a constraint, not an estimate.
Together with the Python notebook this is the full "from-scratch and package" treatment, and the two halves do different jobs: the hand-built binomial weights (validated to machine precision) show how the operator works, and the package supplies the estimator that shows where to stop — agreeing with the from-scratch GPH implementation in the Python companion to three decimals.
Fractional integration is the same long-memory mathematics as ARFIMA in the volatility-persistence arc; here it serves ML as a stationarity-preserving-memory feature transform. This is the package mirror of fracdiff_python.ipynb.