Replicating Functional PCA for Crude Oil Futures Curves

# VAR FPCA analysis paper replication


# VAR FPCA analysis paper replication












I've been trying to replicate the following publication: CONSISTENT FUNCTIONAL PCA FOR FINANCIAL TIME-SERIES, Sebastian Jaimungal, Eddie K. H. Ng, 2007

but I havent been able to get the same results estimating the $\beta_{k}$ parameters.

First, I got the data from Bloomberg (CL Comdty contracts) from the period specified in the paper but I don't get how is the $\tau$ parameter specified. I using it as number of years to maturity but with the 70 contracts available in bbg seems like I won't fit the $\tau$ that appears in the publication's charts.

Also, assuming I'm correct, when i fit the $\beta_{k}$ using MLS with the common formula ($\beta_{k} = (X^{T}X)^{-1}Xy)$ I get, for the first beta, a similar chart but with a different scale:

I guess that the parameter calculation is correct, because when i plot the curve against the data it seems very close and the error is quite small, but cant get the same $ \beta_{k}$ scale:

could it be maybe that the betas are calibrated with the specifict tenors of $\phi$ function?

code:

```
# Libraries:

        install.packages("vars")  #If not already installed
        install.packages("astsa") #If not already installed
        install.packages("fpca")  #If not already installed
        install.packages("matrixcalc",  repos="http://R-Forge.R-project.org")  #If not already installed

    library(vars)
    library(astsa)
    library("XLConnect")

    rm(list=ls())
    "%^%" <- function(S, power)
    { 
       with(eigen(S), vectors %*% (values^power * t(vectors))) 
    }

    data <-readWorksheetFromFile("NYMEX.xlsx",sheet=1,startRow=1,endRow=1000,startCol=2,colTypes = 'numeric')
    tau <-readWorksheetFromFile("NYMEX.xlsx",sheet=2,startRow=1,endRow=1000,startCol=2,colTypes = 'numeric')

#Parameters

    # A_k parametros
    ak <- c(4,2,1,0.2)
    #Beta Matrix
    beta_hat <- matrix(data =NA, nrow = length(data[,1]), ncol = 5)
    # Matrix Price
    price_hat <- matrix(data = NA, nrow = length(data[,1]), ncol = length(data[1,]))

#Beta

for(i in 1:length(data[,1]))
{
#Clean from NA 
    count = 0
    for(j in seq(1,length(data[1,]),by=2))
    {
        if( (is.na(data[i,j]) == FALSE) & (as.numeric(tau[i,j]) > 0))
        {
        count = count + 1 
        }
    }

#Vectors Length
    tau_tmp <- c(1:count)*0
    price_tmp <- c(1:count)*0

    count = 0
    for(j in seq(1,length(data[1,]),by=2))
    {
        if( (is.na(data[i,j]) == FALSE) & (as.numeric(tau[i,j]) > 0))
        {
        count = count +1
        tau_tmp[count] <- as.numeric(tau[i,j])/360 
        price_tmp[count] <- log(as.numeric(data[i,j]),base =exp(1))
        #price_tmp[count] <- as.numeric(data[i,j])
        }
    }

#Phi Length
    phi_tmp <- matrix(nrow =5, ncol = length(tau_tmp))
    #price_tmp <- (price_tmp)/sd(price_tmp)

#Phi Values
    phi_tmp[1,] <- 1
    for(j in 1:length(tau_tmp))
    {
        for (k  in 1: length(ak))
        {   
            phi_tmp[k+1,j] <- (1- exp(-ak[k]*tau_tmp[j]))/ak[k]
        }

    }
#Beta (beta_hat = (x'*x)^-1*x'*y)
    A <- ginv(t(phi_tmp) %*% phi_tmp) %*% t(phi_tmp)
    beta_hat[i,] <- t(price_tmp%*%A)

    for(j in 1:length(phi_tmp[1,]))
    {
        price_hat[i,j] <- sum(beta_hat[i,]*phi_tmp[,j])
    }

    }

#Plot beta
colnames(beta_hat) <- c("Beta1","Beta2","Beta3","Beta4","Beta5")
plot.ts(beta_hat,base = exp(1) , main = "Betas", xlab = "")

#VAR(1)
fitvar2 = VAR(beta_hat, p=1, type="both")

#Plot E
e_hat <- residuals(fitvar2)
A <- e_hat[,1]
#Lag E
B <- c(1:length(e_hat[,1]))
for (i in 1:length(e_hat[,1]))
{
    if (i == 1)
    {
        B[i] <- NA
    }
    if (i >1)
    {
        B[i] <- e_hat[i-1,1]
    }       
}
#Model against data
b <- 857
ph_plot <- na.omit(price_hat[b,])
pr_plot <- na.omit(log(as.numeric(data[b,]),base=exp(1)))
tau_plot <- na.omit(as.numeric(tau[b,])/360)
plot(ph_plot,col='red')
lines(pr_plot,col='blue')

#FPCA
R <- 0.05
space <- 0.05
#Parameters
phi_tmp <- matrix(data = 0, nrow =5, ncol = length(seq(R,7, by=space)))
N= length(e_hat[,1])
mat = seq(R,7, by=space)
    phi_tmp[1,] <- 1
    for(j in 1:length(seq(R,7, by=space)))
    {
        for (k  in 1: length(ak))
        {   
            phi_tmp[k+1,j] <- (1- exp(-ak[k]*mat[j]))/ak[k]
            #phi_tmp[k+1,j] <- (1- exp(-ak[k]*j))/ak[k]

        }
    }

W = phi_tmp %*% t(phi_tmp)
uk = eigen((1/N)*(W %^% 0.5) %*% t(e_hat) %*% e_hat %*% (W %^% 0.5))
z <- (W %^% -0.5) %*% uk$vectors

epsilon <- matrix( data = 0, nrow = 5, ncol = length(phi_tmp[1,]))

for (i in 1:length(phi_tmp[1,]))
{
    epsilon[,i] = t(z) %*% phi_tmp[,i]
}

plot(mat,epsilon[1,], ylim = c(-1,1))
lines(mat,epsilon[2,])
lines(mat,epsilon[3,])

print(uk$values/sum(uk$values))
yt = ginv(z) %*% t(beta_hat)
plot.ts(t(yt))
```