#################################################################### ## ## Script for the Oil Shock model originally developped by WebHubbleTelescope ## see the blog http://www.mobjectivist.com for more details ## ## Adapted by Khebab (April, 2007): ## http://www.theoildrum.com/node/2376 ## ## version 1.0 ## #################################################################### rm(list=ls(all=TRUE)) # clean up graphics.off() # closing graphs ## ## Load matlab emulator package and gremisc (to install) ## in the menu bar: packages/install pacakage(s) ## library(utils) library(matlab) ################################################################### ## ## Safe array retrieval function ## ################################################################### Get <- function(Arr, I) { if (I > 0 && I <= length(Arr)) { return(Arr[I]); } else { return(0.0); } } ################################################################### ## ## Discretized convolution function ## ################################################################### Convolve <- function(A, B) { Total <- length(A) + length(B) + 1; C <- matrix(0, nrow= Total, ncol= 1); for (J in 1:Total) { V <- 0.0; for (I in 1:J) { V <- V + Get(A, I) * Get(B, J - I); } C[J] <- V; } return(C); } ################################################################### ## ## Discretized exponential function ## ################################################################### Exponential <- function(L, Alpha) { R <- matrix(0, nrow= L, ncol= 1); for(I in 1:L) { R[I] <- Alpha * exp(-Alpha * (I - 1)); } R <- R / sum(R); return(R); } ################################################################### ## ## Fills in data between discrete years ## ################################################################### Expander <- function(Arr, Expansion_Factor) { R <- matrix(0, nrow= length(Arr) * Expansion_Factor, ncol= 1); for(I in 1:(length(Arr) * Expansion_Factor)) { R[I] <- Arr[ceil(I / Expansion_Factor)]; } return(R); } ################################################################### ## ## Safe function for reading forcing function ## ################################################################### Discovery_Window <- function(Discovery, Year, Expansion_Factor) { Start <- Discovery$year[1]; Y <- Year - Start; if (Y < 0.0) { return(0.0); } f <- Discovery$value[floor(Y * Expansion_Factor)+1]; return(f); } ################################################################### ## ## Safe function for reading forcing function ## ################################################################### Production_Window <- function(Production, Year) { Start <- Production$year[1]; Y <- floor(Year - Start); if (Y < 0.0) { return(NA); } if (Y > length(Production$year) - 1) { return(NA); } f <- Production$value[Y +1]; return(f); } ################################################################### ## ## Generates the mature discovery by repeated convolution ## ################################################################### Phasing <- function(Strikes, Phases, Time_Span, Expansion_Factor) { Avg_Fallow <- Exponential(Time_Span, Phases[1] / Expansion_Factor); Avg_Build <- Exponential(Time_Span, Phases[2] / Expansion_Factor); Avg_Mature <- Exponential(Time_Span, Phases[3] / Expansion_Factor); startYear <- 1; value <- Expander(Strikes$value, Expansion_Factor); win.graph(); plot(1:100,value[1:100], xlab="", ylab="Gb/year"); value <- Convolve(value, Avg_Fallow); lines(1:100,value[1:100],col="red"); valueB <- Convolve(value, Avg_Build); lines(1:100,valueB[1:100],col="blue"); value <- Convolve(valueB, Avg_Mature); fSum= sum(value); lines(1:100,value[1:100],col="purple"); year <- seq(0,length(value)-1) / Expansion_Factor + Strikes$year[1]; growth <- matrix(0,nrow= length(value), ncol= 1); idx <- find(year >= 1996 & year <= 2025); title(sprintf("Oil Discoveries: %6.0f Gb",fSum)); idx <- find(year <= Strikes$year[1] + Time_Span); year <- year[idx]; value <- value[idx]; growth <- growth[idx]; Discovery <-data.frame(year ,value, growth); return(Discovery); } ################################################################### ## ## Compute ## ################################################################### Compute <- function(Discovery, Production, Time_Span, Expansion_Factor) { C <- 0; # -- Reserve DT <- 1; # -- time step Start <- Discovery$year[1]; Finish <- Time_Span + Start; V <- 0.0; #-- Volume of oil extracted Data <- 0.0; n <- 1; Time <- Start; NParticle <- 200; fSigma <- 0.25*0.25; fKMax <- 1/20; CParticle <- matrix(0, nrow= NParticle , ncol= 1); PParticle <- matrix(0, nrow= NParticle , ncol= 1); ParticleWeight <- matrix(0, nrow= NParticle , ncol= 1); CParticle2 <- matrix(0, nrow= NParticle , ncol= 1); PParticle2 <- matrix(0, nrow= NParticle , ncol= 1); Particles2 <- matrix(0, nrow= NParticle , ncol= 1); K <- matrix(0,nrow= 1,ncol= 1); R <- matrix(0,nrow= 1,ncol= 1); D <- matrix(0,nrow= 1,ncol= 1); win.graph(); par(lwd= 2); dev.set(4); while(Time <= Finish) { #-- Integration of discovery window against extraction rate fDiscovery <- Discovery_Window(Discovery, Time, Expansion_Factor); fObservedProduction <- Production_Window(Production, Time); if (n == 1) { fMeanRate <- 0.05; Particles <- rnorm(NParticle, mean= fMeanRate, sd= 0.03) ; } else { Particles <- Particles + rnorm(NParticle, mean= 0, sd= 0.03) ; } if (!is.na(fObservedProduction)) { for(m in 1:NParticle) { Rate <- Particles[m]; CParticle[m] <- CParticle[m] + fDiscovery * DT - CParticle[m] * Rate * DT; PParticle[m] <- Rate * CParticle[m]; ParticleWeight[m] <- exp(-0.5 * (fObservedProduction - PParticle[m]) * (fObservedProduction - PParticle[m])/fSigma) } ### Weight ### fSum <- sum(ParticleWeight); for(m in 1:NParticle) { ParticleWeight[m]<- ParticleWeight[m] / fSum; } ### Importance Resampling ### vNewIndex <- sample(c(1:NParticle), NParticle, replace = TRUE, prob = ParticleWeight); for(m in 1:NParticle) { Particles2[m]<- Particles[vNewIndex[m]]; CParticle2[m]<- CParticle[vNewIndex[m]]; PParticle2[m]<- PParticle[vNewIndex[m]]; } Particles <- Particles2; CParticle <- CParticle2; PParticle <- PParticle2; fMeanRate <- mean(Particles); K[n] <- fMeanRate; } else { ### Prediction ### for(m in 1:NParticle) { CParticle[m] <- CParticle[m] + fDiscovery * DT - CParticle[m] * fMeanRate * DT; PParticle[m] <- fMeanRate * CParticle[m]; } } value[n] <- mean(PParticle); year[n] <- Time; K[n] <- fMeanRate; R[n] <- mean(CParticle); D[n] <- fDiscovery; if (year[n] > 1970 & year[n] < 1980) { idx <- find(year == 1970); fKMax= max(K[idx:n]); } plot(year[1:n],value[1:n], panel.first=grid(8,8), type="l", xlab="", ylab="Gb/year", xlim=c(1930,2040), ylim=c(0,50)); lines(year[1:n],D[1:n], col="blue"); text(x= year[which.max(value)], y= max(value)+0.5, sprintf("PO= %6.3f", year[which.max(value)]), col= 1, cex= .8) title(sprintf("E(%d)= %6.3f %% (%6.2f Gb)", Time,fMeanRate*100,value[n])); n <- n + 1; Time <- Time + DT; } Results <-data.frame(year ,value, R, D, K); return(Results); } ################################################################### ## ## Main ## #################################################################### Phases <- c(1/3, 1/3, 1/3); # 1/lambda # Computation over 200 years Time_Span <- 200; # Calculation precision Expansion_Factor <- 1; # # World oil production since 1900 up to 2006 in Gb/year (crude oil + condensate) # 1900-1959: API Facts and Figures Centennial edition 1959. # 1960-2005: http://www.eia.doe.gov/emeu/international/oilproduction.html # temp <- scan() 2.0000000e-001 1.6740000e-001 1.8180000e-001 1.9490000e-001 2.1800000e-001 2.1510000e-001 2.1330000e-001 2.6420000e-001 2.8560000e-001 2.9870000e-001 3.2780000e-001 3.4440000e-001 3.5240000e-001 3.8530000e-001 4.0750000e-001 4.3200000e-001 4.5750000e-001 5.0290000e-001 5.0350000e-001 5.5590000e-001 6.8890000e-001 7.6600000e-001 8.5890000e-001 1.0157000e+000 1.0143000e+000 1.0679000e+000 1.0968000e+000 1.2630000e+000 1.3250000e+000 1.4860000e+000 1.4100000e+000 1.3730000e+000 1.3100000e+000 1.4420000e+000 1.5220000e+000 1.6540000e+000 1.7920000e+000 2.0390000e+000 1.9880000e+000 2.0860000e+000 2.1500000e+000 2.2210000e+000 2.0930000e+000 2.2570000e+000 2.5930000e+000 2.5950000e+000 2.7450000e+000 3.0220000e+000 3.4330000e+000 3.4040000e+000 3.8030000e+000 4.2830000e+000 4.5190000e+000 4.7980000e+000 5.0180000e+000 5.6260000e+000 6.1250000e+000 6.4390000e+000 6.6080000e+000 7.1340000e+000 7.6614E+00 8.1943E+00 8.8878E+00 9.5375E+00 1.0286E+01 1.1070E+01 1.2030E+01 1.2917E+01 1.4100E+01 1.5221E+01 1.6750E+01 1.7710E+01 1.8666E+01 2.0323E+01 2.0338E+01 1.9283E+01 2.0929E+01 2.1794E+01 2.1958E+01 2.2875E+01 2.1754E+01 2.0469E+01 1.9520E+01 1.9440E+01 1.9889E+01 1.9703E+01 2.0524E+01 2.0685E+01 2.1440E+01 2.1849E+01 2.2108E+01 2.1977E+01 2.1977E+01 2.1988E+01 2.2261E+01 2.2750E+01 2.3251E+01 2.3977E+01 2.4426E+01 2.4035E+01 2.4955E+01 2.4813E+01 2.4444E+01 2.5269E+01 2.6360E+01 2.6857E+01 2.68129E+01 value <- matrix(temp, ncol=1 , byrow= F); year <- 1900+c(1:length(value))-1; vWorldProduction <-data.frame(year ,value); #### Nominal oil prices (constant 2000 dollars) ### value <- scan() 13.5500 8.1200 12.1100 9.8300 14.1900 13.4300 14.9300 15.6000 15.2400 13.9500 13.8200 14.7000 13.8600 13.1800 13.0700 11.0800 10.9000 16.1500 15.6900 14.1800 13.4900 12.5000 12.2200 13.6900 13.6300 13.6800 13.4800 12.8100 13.6600 13.5500 12.1800 11.4200 11.2900 11.1500 11.0100 10.8300 10.5000 10.2300 9.8200 9.3200 8.7900 10.5000 11.2600 14.0500 44.5500 40.6600 41.2600 41.6200 39.5500 78.4600 82.1500 71.4700 62.3600 54.7300 51.1200 48.4600 24.7900 30.6400 23.9000 27.7500 34.4400 27.8400 26.0900 22.3300 20.3800 21.3100 25.0800 22.7400 15.2000 20.7100 31.8000 26.4400 26.4700 29.6100 38.2700 50.0000 year <- 1930+c(1:length(value))-1; vPrices <-data.frame(year ,value); #### Proven Reserves (BP, 2006) ### value <- scan() 6.671E+02 6.875E+02 7.168E+02 7.274E+02 7.616E+02 7.704E+02 8.773E+02 9.089E+02 9.962E+02 1.004E+03 1.001E+03 1.006E+03 1.011E+03 1.012E+03 1.018E+03 1.027E+03 1.049E+03 1.059E+03 1.070E+03 1.090E+03 1.115E+03 1.140E+03 1.173E+03 1.188E+03 1.194E+03 1.201E+03 #### Corrected Proven Reserves (Middle-East reserve jumps removed) ### valueC <- scan() 667.1000 687.5000 698.5000 709.1000 720.1000 728.9000 739.9000 750.9000 761.9000 769.7000 766.7000 771.7000 776.7000 777.7000 783.7000 792.7000 814.7000 824.7000 835.7000 855.7000 866.7000 877.7000 888.7000 903.7000 909.7000 916.7000 year <- c(1980:2005); ProvenReserves <-data.frame (year ,value, valueC); #-- Data transcribed from ASPO discovery curve #-- http://wolf.readinglitho.co.uk/mainpages/discoveries.html #-- GBBls/year starting year 1932 # discovery at 2Gb/year between 1912 and 1931 (~average for the US) # value <- scan() 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.0 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.0 1.5 9.0 4.0 4.5 3.0 5.0 5.5 42.0 42.5 52.0 18.0 16.0 5.0 3.0 7.4 7.6 7.0 49.0 53.0 54.0 19.5 16.0 26.0 20.0 28.0 24.0 35.0 40.0 40.5 42.0 44.0 50.0 41.0 50.0 48.5 47.5 32.0 30.5 30.4 29.5 40.5 37.0 38.5 24.0 26.5 31.0 37.0 38.0 36.0 26.0 24.0 20.0 20.0 22.0 21.0 20.5 18.0 16.2 17.5 16.5 20.0 17.0 16.0 9.0 8.5 9.0 9.0 9.5 14.0 18.0 17.5 13.0 9.0 9.0 # from 2005 to 2031, we discover 11Gb/year value[93:120]= 11; year <- c(1912:2031); ASPODiscovery <-data.frame(year ,value); #### Discovery Logistic Decay up to 2033 ### valueHubbert <- scan() 8.815E+00 8.428E+00 7.989E+00 7.643E+00 7.245E+00 6.893E+00 6.539E+00 6.232E+00 5.924E+00 5.613E+00 5.352E+00 5.038E+00 4.827E+00 4.562E+00 4.295E+00 4.081E+00 3.866E+00 3.649E+00 3.486E+00 3.269E+00 3.104E+00 2.940E+00 2.775E+00 2.664E+00 2.498E+00 2.387E+00 2.220E+00 2.108E+00 1.997E+00 1.885E+00 value[93:122]= valueHubbert; year <- c(1912:2033); ASPODiscovery <-data.frame(year ,value); plot(ASPODiscovery$year, ASPODiscovery$value, xlab="Year", ylab="Gb/year", type="h"); title('ASPO Discoveries'); Discoveries <- Phasing(ASPODiscovery, Phases, Time_Span, Expansion_Factor); dev.set(2); lines(Discoveries$year, Discoveries$value, xlim=c(1930,2040), col="blue") Results <- Compute(Discoveries, vWorldProduction, Time_Span, Expansion_Factor); dev.set(4); lines(vWorldProduction$year, vWorldProduction$value, type="o", col="red"); dev.set(2); lines(vWorldProduction$year, vWorldProduction$value, col="red"); lines(Results$year, Results$value); write.csv(Results, file = "Results.csv") win.graph(); idx= find(vWorldProduction$year >= 1980 & vWorldProduction$year <= 2005); plot(Results$year, Results$K^-1, type="l",col="blue", xlab="", ylab="R/P Ratio (years)", xlim=c(1930,2020), ylim=c(0,80)); lines(ProvenReserves$year, ProvenReserves$value/vWorldProduction$value[idx], type="p"); lines(ProvenReserves$year, ProvenReserves$valueC/vWorldProduction$value[idx], type="p", col= "red"); lines(vPrices$year, vPrices$value);