Created
December 27, 2011 01:14
-
-
Save robbrit/1522394 to your computer and use it in GitHub Desktop.
Canadian Oil Production
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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