Given below is a brief outline of how you might go about modelling long run oil prices based on the work by Pindyck (1999) on "The long-run evolution of energy prices", published in the 'Energy Journal'. The data for the example are taken from the BP Statistical Review of Energy, which is freely available online. This includes price data since 1861, deflated to 2010 prices, and are a little different from those used by Pindyck (1999). The code writes the data into 'R' (as if typing in by hand - see other examples for alternative data entry methods) then forms a state space object to be estimated by maximum likelihood. Specifically, prices are allowed to follow a first order autoregressive process but with additional controls for time that are partially stochastic and partially deterministic:

Pt=ρPt-1+ β12t +β3t21t+ φ2tt +εt(1)

Where α are the prices in period 't' and t is the time index. α, β1, β2 and β3 are static parameters. φ1t and φ2t are time varying parameters to be estimated using the same state equations as Pindyck (1999):

φ1t= c1φ1t-1t(2)

φ2t= c2φ2t-1t(2)

c1 and c2 are static parameters to be estimated. The code below describes how to define and estimate this system using 'R', but noting that the final estimated equations are a little different from those stated above - for further details refer to Pindyck (1999).

Click the [+] to expand or collapse a section.


[+] Data entry/preperation

    Step 1: Enter the data as a variable named 'oil'

    oil<-c(11.84371818168, 22.84145649324, 55.5602995781513, 111.916194124933, 
    93.4939119816469, 55.47210862644, 37.4475579242914, 59.224633621758, 
    59.387786882424, 66.2917459127115, 78.67612792116, 65.98642986936, 
    33.17449633542, 22.4575664681435, 26.6978062908, 52.20904341312, 
    49.353861351465, 26.7796386472469, 20.0445434532514, 21.3787031217517, 
    19.353352299691, 17.5530404578593, 23.3076086665714, 20.30351688288, 
    21.27035102016, 17.16130593672, 16.19447179944, 21.27035102016, 
    22.72060222608, 21.02864248584, 16.19447179944, 13.53567792192, 
    15.46934619648, 21.0844213783754, 35.5021495209216, 30.8033356137408, 
    20.6225721481824, 23.7551147529696, 33.6748330014624, 31.0643808308064, 
    25.0603408382976, 20.0804013127385, 22.72060222608, 20.78693395152, 
    14.98592912784, 17.6447230053599, 16.7814782399314, 17.40301447104, 
    16.9195974024, 14.2176412866086, 14.2176412866086, 16.6528845369434, 
    20.874827964, 17.5600932912, 13.7364813312, 21.952133632, 26.5092423648, 
    28.646670888, 25.3289270424, 33.3979584888, 21.068679748, 20.9308029392, 
    17.1138634928, 18.225883676, 20.8847059008, 23.1496099616, 16.31385964, 
    14.8814715816, 16.1533922296, 15.5328702704, 9.302159932, 13.8809870424, 
    11.2670407424, 16.27569984, 15.4003576392, 17.1439334152, 17.9110762288, 
    17.4773846504, 16.002911784, 15.8494439712, 16.8738710592, 15.9117207648, 
    15.121747488, 14.986552768, 12.711792576, 12.491992128, 18.536068336, 
    18.0124940744, 16.2746967824, 15.4780727976, 14.3408017296, 14.038772364, 
    15.718675788, 15.6429231336, 15.7060503456, 15.4787923824, 14.703734136, 
    15.6794475136, 15.5433805696, 13.974554872, 13.10952672, 12.968226432, 
    12.81515112, 12.63852576, 12.426575328, 12.08902464, 11.747548944, 
    11.276547984, 10.699571808, 10.102970592, 12.0450645504, 12.9300230208, 
    16.1487694424, 51.2340452592, 46.7387117912, 49.039922176, 50.0527486848, 
    46.8887288343558, 94.9414622589532, 97.4636223300971, 86.1908919691969, 
    74.500583626943, 64.6943253012048, 60.4008823869105, 55.8515182156134, 
    28.7093802919708, 35.3861879525211, 27.5083112910837, 32.0509142076194, 
    39.5834542565937, 32.0253903316535, 30.0286838606529, 25.6108419869564, 
    23.273012707498, 24.3476975634068, 28.7245883145341, 25.939272647352, 
    17.0105906996564, 23.5203078316735, 36.0836450101928, 30.0967652044133, 
    30.3305784698432, 34.1668901911609, 44.1710579142403, 60.8737874654992, 
    70.4615758556548, 76.1296451253986, 98.4995490303284, 62.6828440652605, 
    79.4955336)
    	

    Step 2: Organise the data into a time series object, take natural logs then define a restricted time window
    oil<-ts(oil,start=1861)
    
    oil<-log(oil)
    
    oil_pindyck<-window(oil, start=1870,end=1996)
    
    plot(oil_pindyck, main="Log of oil prices in constant 2010 dollars.")
    	

[+] Setting up the model

    Step 1: Generate the remaining estimation variables

    P<-oil_pindyck[2:length(oil_pindyck)]
    P_lag<-oil_pindyck[1:(length(oil_pindyck)-1)]
    
    t<-1:length(P)
    t_squared<-t*t
    
    	

    Step 2: Load the 'tsa3.rda' data archive containing the Kalman filter code. Note that you will need this file to be in the current working directory, you may need to change the working directory. The file can be found on the website for the book "Time series analysis and its applications" by Robert Shumway and David Stoffer.
    load("tsa3.rda")
    	

    Step 3: Write the state space system and estimate
    y<-P
    num=length(y)
    A<-array(1, dim=c(1,2,num))
    A[,2,]<-t
    
    mu0 = rep(0.1,2)
    Sigma0=diag(100,2)
    Phi=diag(0,2)
    
    # Function to Calculate Likelihood 
    Linn=function(para){
    
    cR=para[1]
    
    cQ<-diag(0,2)
    cQ[1,1]=para[2]
    cQ[2,2]=para[3]
    
    input<-cbind(1,P_lag)
    
    Gam<-matrix(1,1,2)
    Gam[1,1]<-para[4]
    Gam[1,2]<-para[5]
    
    Phi[1,1]<-para[6]
    Phi[2,2]<-para[7]
    
      kf = Kfilter1(num,y,A,mu0,Sigma0,Phi,0,Gam,cQ,cR,input)  
     return(kf$like) 
      list(kf=kf)
     }
    (est = optim(rep(-0.1,7), Linn, NULL, method="BFGS", hessian=FALSE, 
    control=list(trace=1,REPORT=1,maxit=5000,reltol=sqrt(.Machine$double.eps) ))) 
    

[+] Extracting the main results

    Step 1: Extract the estimated parametrs and run the Kalman filter at the estimates

    para<-est$par
    
    cR=para[1]
    
    cQ<-diag(0,2)
    cQ[1,1]=para[2];cQ[2,2]=para[3]
    
    input<-cbind(1,P_lag)
    
    Gam<-matrix(1,1,2)
    Gam[1,1]<-para[4];Gam[1,2]<-para[5]
    
    Phi[1,1]<-para[6];Phi[2,2]<-para[7]
    
      kf = Kfilter1(num,y,A,mu0,Sigma0,Phi,0,Gam,cQ,cR,input)
    

    Step 2: Plot the stochastic/time varting terms
    layout(rbind(1,2))
    plot(ts(kf$xf[1,,], start=1871), main="Phi 1")
    plot(ts(kf$xf[2,,], start=1871), main="Phi 2")
    

[+] A simple 20 year forecast

    Step 1: Store main parameters and create forecast variables

    c1=para[4]
    c2=para[5]
    c3=para[6]
    c4=para[7]
    P_lag_forecast<-c(P_lag,rep(NA,20))
    t_forecast<-c(t,rep(NA,20))
    phi1_forecast<-c(kf$xf[1,,],rep(NA,20))
    phi2_forecast<-c(kf$xf[2,,],rep(NA,20))
    P_hat<-c1+(c2*P_lag_forecast)+phi1_forecast+(phi2_forecast*t_forecast)
    	

    Step 2: Recursively update the forecast variables
    for (i in (length(P_hat)-20):(length(P_hat)) )
    {
    P_lag_forecast[i] <- P_hat[i-1]
    t_forecast[i] <- (t_forecast[i-1]+1)
    phi1_forecast[i]<-(c3*phi1_forecast[i-1])
    phi2_forecast[i]<-((c4*phi2_forecast[i-1]) )
    
    P_hat[i]<-c1+(c2*P_lag_forecast[i])+(phi1_forecast[i])+(phi2_forecast[i]*t_forecast[i])
    }
    	

    Step 3: Convert the forecast into a time series object, then plot the forecast along with the actual data (where available)
    P_hat<-ts(P_hat, start=1871)  # noting that starts in 1871 due to one period lag
    plot(P_hat, main="Comparison of forecast with hold-out sample",
    xlab="(Solid line=forecast; Dashed line=actual.)")
    abline(v=1996)
    lines(oil,lty=2)