Create a gist now

Instantly share code, notes, and snippets.

What would you like to do?
Canadian Oil Production
oil = read.csv("oil.csv", header = T)
# convert cubic metres to barrels
total = ts(oil$Total * 6.28981077 / 1000, start = 1971, end = 2010)
t = 1:length(total)
linear = lm(total ~ t)
expn = lm(log(total) ~ t)
# do a linear model
png("oil.png")
plot(total, xlab = "Year", ylab = "Production Mbbl")
title("Linear Model")
# need to offset intercept so that we can see the line on the plot
abline(linear$coefficients[1] - linear$coefficients[2] * 1971, linear$coefficients[2], col = "blue")
dev.off()
print(linear$coefficients[2])
# do an exponential model
png("oilx.png")
plot(as.vector(total), xlab = "Year from 1971", ylab = "Production Mbbl")
summary(expn)
# subtract 1 from t so that we start at 0
expn_plot = total[1] * exp(expn$coefficients[2] * (t - 1))
lines(expn_plot, col = "blue")
title("Exponential Model")
print(expn$coefficients[2])
# do a linear and exponential model with 1970's bubble excluded
total2 = ts(total[11:length(total)], start = 1981, end = 2010)
t2 = 1:length(total2)
linear2 = lm(as.vector(total2) ~ t2)
expn2 = lm(log(total2) ~ t2)
summary(linear2)
png("oilt.png")
plot(as.vector(total2), xlab = "Year from 1981", ylab = "Production Mbbl")
abline(linear2$coefficients[1], linear2$coefficients[2], col = "blue")
title("Linear from 1981")
dev.off()
png("oilxt.png")
plot(as.vector(total2), xlab = "Year from 1981", ylab = "Production Mbbl")
expn_plot2 = total2[1] * exp(expn2$coefficients[2] * (t2 - 1))
lines(expn_plot2, col = "blue")
title("Exponential from 1981")
dev.off()
summary(expn2)
print(expn2$coefficients[2])
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment