Here is the result of my recent workouts on the returns data you provided me with and that I performed as an exercise in relation with my long journey through the vast, dry and dull R territory. Having jumped that hurdle I now can make out an oasis in the far end of the desert; a mirage for sure.
So here is the thing. You'll find bellow stats related to grains. I start with a brief summary of each individual distributions for the full period (1992-2011). The norm plot, box plot, density and QQ plot provide a convenient overview of the distribution at a glance. I then compute the 4 first moments with corresponding standard error and 95% confidence interval. Finally, I take a close, longitudinal look with the time series of the 52-week rolling \( x^{th} \) moment over the whole period. See my comments in the code.
options(digits = 4, width = 70)
library(boot)
library(tseries)
library(zoo)
library(PerformanceAnalytics)
library(multicore)
returns.df = read.csv(file = "~/Documents/annual/returnsCME.csv", sep = ";",
header = TRUE, stringsAsFactors = FALSE)
returns.df[is.na(returns.df)] = 0
returns.z = zooreg(data = returns.df[, -1], start = c(1992, 41), end = c(2011,
12), frequency = 52)
grains.z = returns.z[, c("corn", "oats", "wheat", "rice", "soybeans", "soymeal",
"soyoil")]
wget.and.source <- function(url) {
fname <- tempfile()
download.file(url, fname, method = "wget")
source(fname)
unlink(fname)
}
wget.and.source("https://gist.github.com/bautheac/b2e8269868a881b66a69/raw/")
rollMomentsGrains.ls = rollingMoments.fun(grains.z, startYear = 1992, startWeek = 41,
width = 52)
mu95CI(grains.z, "corn")
## muhat standard error lower bound upper bound
## [1,] -0.001183 0.001214 -0.00361 0.001245
The standard error of the mean can be computed analytically with no error. \( \widehat{\mu} \) = \( \widehat{\sigma} /\sqrt{t} \) So far so good.
sd95CI(grains.z, "corn")
## sdhat standard error lower bound upper bound
## [1,] 0.03768 0.001301 0.03508 0.04028
Here it got tedious. The analytic solution for the standard error of sigma is only an approximation ( \( \widehat{\sigma} \) \( \approx \widehat{\sigma}/\sqrt{2\times t} \) ). For increased precision, I use the bootstrapping method. The times series is resampled 999 times; for each sample, sigma is computed. The standard error is the standard deviation of the 999 sigmas estimates.
I proceed similarly for skewness and kurtosis.
skew95CI(returns.df, "corn")
## skewhat standard error lower bound upper bound
## [1,] -0.1425 0.2119 -0.5663 0.2813
kurt95CI(returns.df, "corn")
## kurthat standard error lower bound upper bound
## [1,] 2.726 0.5841 1.557 3.894
Here it stepped up from tedious to daunting. No issue with the mean still, I used a convient “rollApply” R function that allows to apply a particular function to a times series in a rolling fashion over a given period according to a specified width (52 here). I only had to write a function that returns the mean as well as the 95% CI drawn from the computed standard error. However, sigma and higher moments call for bootstrap which implied 999 operations for a single 52-week period; repeated over the whole period (1993-2011) for the 3 moments and the 28 variables, it boils down to roughly 90 million computations. In other words, blaze of glory for the old home computer I'm working with here. I first tried a decent CPU overcloking and rapidly reached the stability limits of the system. I then found out about this R packages that allows to have all the CPU cores work simultaneously for heavy computational effort; ray of hope. I got down in the office where there is a more recent computer with 4 CPU cores and found out the damn R package doesn't run on windows but on linux only. Got linux installed along side windows in dual boot; second ray of hope, vanished early. No more convenient rolling function with this gofast thing and I had to rewrite the all damn code. After 2 days searching the web for a convenient solution, insulting the monitor and injecting coffee I ended with a loop that creates a list of ranges with the first range containing observations 1 to 52, the second containing observations 2 to 53 and so on; I than apply the fast routine to every single range. It still requires half an hour to compile but it's not the all day and the output offers a decent level of precision. There is a resilient issue however related to the data. the best code I might write doesn't prevent “garbage in, garbage out”: the returns time series show a lot of seros and NAs. Glancing at the plots you'll see some of them only display data after 2000 and even after 2007 for some ernegy commodities time series. That's because I wrote the plotting functions such that they find the first row where all entries are non zeros and proceeds from there.
head(rollMomentsGrains.ls[["corn"]][["mu"]])
## muhat.lower muhat muhat.upper
## 1993(41) -0.043488 -1.000e-02 0.02349
## 1993(42) -0.033373 -3.464e-18 0.03337
## 1993(43) -0.004394 3.000e-02 0.06439
## 1993(44) 0.015187 5.000e-02 0.08481
## 1993(45) 0.034791 7.000e-02 0.10521
## 1993(46) 0.054416 9.000e-02 0.12558
plotRollingMoments(data = rollMomentsGrains.ls, commo = "corn", statistic = "mu")
head(rollMomentsGrains.ls[["corn"]][["sd"]])
## sdhat.lower sdhat sdhat.upper
## 1993(41) 0.1161 0.1207 0.1254
## 1993(42) 0.1159 0.1203 0.1248
## 1993(43) 0.1197 0.1240 0.1283
## 1993(44) 0.1212 0.1255 0.1299
## 1993(45) 0.1229 0.1269 0.1310
## 1993(46) 0.1243 0.1283 0.1323
plotRollingMoments(data = rollMomentsGrains.ls, commo = "corn", statistic = "sd")
head(rollMomentsGrains.ls[["corn"]][["skew"]])
## skewhat.lower skewhat skewhat.upper
## 1993(41) 0.158439 1.0682 1.978
## 1993(42) 0.151770 1.0483 1.945
## 1993(43) 0.141053 0.9611 1.781
## 1993(44) 0.083147 0.8874 1.692
## 1993(45) 0.040742 0.8167 1.593
## 1993(46) -0.009616 0.7484 1.507
plotRollingMoments(data = rollMomentsGrains.ls, commo = "corn", statistic = "skew")
head(rollMomentsGrains.ls[["corn"]][["kurt"]])
## kurthat.lower kurthat kurthat.upper
## 1993(41) -0.8845 1.9604 4.805
## 1993(42) -0.8796 1.9772 4.834
## 1993(43) -1.0029 1.4597 3.922
## 1993(44) -1.1131 1.1979 3.509
## 1993(45) -1.1771 0.9644 3.106
## 1993(46) -1.2630 0.7559 2.775
plotRollingMoments(data = rollMomentsGrains.ls, commo = "corn", statistic = "kurt")
mu95CI(grains.z, "oats")
## muhat standard error lower bound upper bound
## [1,] -0.0004046 0.001428 -0.00326 0.002451
sd95CI(grains.z, "oats")
## sdhat standard error lower bound upper bound
## [1,] 0.04432 0.001464 0.0414 0.04725
skew95CI(returns.df, "oats")
## skewhat standard error lower bound upper bound
## [1,] -0.0856 0.1831 -0.4518 0.2806
kurt95CI(returns.df, "oats")
## kurthat standard error lower bound upper bound
## [1,] 1.994 0.5045 0.9852 3.003
head(rollMomentsGrains.ls[["oats"]][["mu"]])
## muhat.lower muhat muhat.upper
## 1993(41) -0.1611 -0.11 -0.05892
## 1993(42) -0.1611 -0.11 -0.05892
## 1993(43) -0.1417 -0.09 -0.03825
## 1993(44) -0.1620 -0.11 -0.05801
## 1993(45) -0.2137 -0.16 -0.10634
## 1993(46) -0.1750 -0.12 -0.06503
plotRollingMoments(data = rollMomentsGrains.ls, commo = "oats", statistic = "mu")
head(rollMomentsGrains.ls[["oats"]][["sd"]])
## sdhat.lower sdhat sdhat.upper
## 1993(41) 0.1772 0.1842 0.1911
## 1993(42) 0.1773 0.1842 0.1911
## 1993(43) 0.1797 0.1866 0.1935
## 1993(44) 0.1808 0.1875 0.1941
## 1993(45) 0.1869 0.1935 0.2000
## 1993(46) 0.1919 0.1982 0.2045
plotRollingMoments(data = rollMomentsGrains.ls, commo = "oats", statistic = "sd")
head(rollMomentsGrains.ls[["oats"]][["skew"]])
## skewhat.lower skewhat skewhat.upper
## 1993(41) -0.050557 0.8303 1.711
## 1993(42) -0.046361 0.8303 1.707
## 1993(43) -0.012446 0.7893 1.591
## 1993(44) 0.005897 0.8161 1.626
## 1993(45) -0.017148 0.7386 1.494
## 1993(46) -0.027926 0.6784 1.385
plotRollingMoments(data = rollMomentsGrains.ls, commo = "oats", statistic = "skew")
head(rollMomentsGrains.ls[["oats"]][["kurt"]])
## kurthat.lower kurthat kurthat.upper
## 1993(41) -0.4188 1.6952 3.809
## 1993(42) -0.4541 1.6952 3.844
## 1993(43) -0.4878 1.4548 3.398
## 1993(44) -0.5171 1.4240 3.365
## 1993(45) -0.5845 1.2024 2.989
## 1993(46) -0.8207 0.8549 2.530
plotRollingMoments(data = rollMomentsGrains.ls, commo = "oats", statistic = "kurt")
mu95CI(grains.z, "wheat")
## muhat standard error lower bound upper bound
## [1,] -0.001463 0.001275 -0.004014 0.001088
sd95CI(grains.z, "wheat")
## sdhat standard error lower bound upper bound
## [1,] 0.0396 0.001169 0.03726 0.04194
skew95CI(returns.df, "wheat")
## skewhat standard error lower bound upper bound
## [1,] 0.1408 0.1708 -0.2008 0.4824
kurt95CI(returns.df, "wheat")
## kurthat standard error lower bound upper bound
## [1,] 1.685 0.4421 0.8005 2.569
head(rollMomentsGrains.ls[["wheat"]][["mu"]])
## muhat.lower muhat muhat.upper
## 1993(41) 0.05861 0.11 0.1614
## 1993(42) 0.09877 0.15 0.2012
## 1993(43) 0.07930 0.13 0.1807
## 1993(44) 0.12914 0.18 0.2309
## 1993(45) 0.10935 0.16 0.2106
## 1993(46) 0.09938 0.15 0.2006
plotRollingMoments(data = rollMomentsGrains.ls, commo = "wheat", statistic = "mu")
head(rollMomentsGrains.ls[["wheat"]][["sd"]])
## sdhat.lower sdhat sdhat.upper
## 1993(41) 0.1804 0.1853 0.1901
## 1993(42) 0.1800 0.1847 0.1894
## 1993(43) 0.1781 0.1828 0.1875
## 1993(44) 0.1785 0.1834 0.1882
## 1993(45) 0.1778 0.1826 0.1875
## 1993(46) 0.1776 0.1825 0.1874
plotRollingMoments(data = rollMomentsGrains.ls, commo = "wheat", statistic = "sd")
head(rollMomentsGrains.ls[["wheat"]][["skew"]])
## skewhat.lower skewhat skewhat.upper
## 1993(41) -0.1070 0.4002 0.9074
## 1993(42) -0.1810 0.3323 0.8456
## 1993(43) -0.1508 0.3650 0.8809
## 1993(44) -0.2481 0.2859 0.8198
## 1993(45) -0.2082 0.3298 0.8679
## 1993(46) -0.1985 0.3530 0.9045
plotRollingMoments(data = rollMomentsGrains.ls, commo = "wheat", statistic = "skew")
head(rollMomentsGrains.ls[["wheat"]][["kurt"]])
## kurthat.lower kurthat kurthat.upper
## 1993(41) -1.073 -0.17293 0.7275
## 1993(42) -1.137 -0.19122 0.7543
## 1993(43) -1.021 -0.07639 0.8682
## 1993(44) -1.032 -0.14940 0.7336
## 1993(45) -1.073 -0.08708 0.8985
## 1993(46) -1.027 -0.06950 0.8883
plotRollingMoments(data = rollMomentsGrains.ls, commo = "wheat", statistic = "kurt")
mu95CI(grains.z, "rice")
## muhat standard error lower bound upper bound
## [1,] -0.001048 0.0009523 -0.002952 0.0008568
sd95CI(grains.z, "rice")
## sdhat standard error lower bound upper bound
## [1,] 0.02957 0.001263 0.02704 0.03209
skew95CI(returns.df, "rice")
## skewhat standard error lower bound upper bound
## [1,] -0.4808 0.3278 -1.136 0.1748
kurt95CI(returns.df, "rice")
## kurthat standard error lower bound upper bound
## [1,] 5.026 1.447 2.133 7.92
Here we have 0s in the time series up to year 2000. It brings some weirdos in the sigma, skewness and kurtosis plots. As a consequence, only refer to post 2001 data in the plots please.
head(rollMomentsGrains.ls[["rice"]][["mu"]])
## muhat.lower muhat muhat.upper
## 1993(41) 0 0 0
## 1993(42) 0 0 0
## 1993(43) 0 0 0
## 1993(44) 0 0 0
## 1993(45) 0 0 0
## 1993(46) 0 0 0
plotRollingMoments(data = rollMomentsGrains.ls, commo = "rice", statistic = "mu")
head(rollMomentsGrains.ls[["rice"]][["sd"]])
## sdhat.lower sdhat sdhat.upper
## 1993(41) 0 0 0
## 1993(42) 0 0 0
## 1993(43) 0 0 0
## 1993(44) 0 0 0
## 1993(45) 0 0 0
## 1993(46) 0 0 0
plotRollingMoments(data = rollMomentsGrains.ls, commo = "rice", statistic = "sd")
tail(rollMomentsGrains.ls[["rice"]][["skew"]])
## skewhat.lower skewhat skewhat.upper
## 2011(12) -0.6636 -0.1850 0.2937
## 2011(13) -0.6280 -0.1635 0.3011
## 2011(14) -0.6533 -0.1936 0.2662
## 2011(15) -0.7469 -0.2798 0.1873
## 2011(16) -0.7340 -0.2792 0.1757
## 2011(17) -0.7825 -0.3308 0.1209
plotRollingMoments(data = rollMomentsGrains.ls, commo = "rice", statistic = "skew")
tail(rollMomentsGrains.ls[["rice"]][["kurt"]])
## kurthat.lower kurthat kurthat.upper
## 2011(12) -1.271 -0.5419 0.187760
## 2011(13) -1.329 -0.6629 0.003213
## 2011(14) -1.271 -0.5606 0.149739
## 2011(15) -1.278 -0.4999 0.277979
## 2011(16) -1.290 -0.5604 0.168999
## 2011(17) -1.349 -0.5939 0.161584
plotRollingMoments(data = rollMomentsGrains.ls, commo = "rice", statistic = "kurt")
mu95CI(grains.z, "soybeans")
## muhat standard error lower bound upper bound
## [1,] 0.001297 0.00107 -0.0008434 0.003437
sd95CI(grains.z, "soybeans")
## sdhat standard error lower bound upper bound
## [1,] 0.03322 0.001062 0.0311 0.03535
skew95CI(returns.df, "soybeans")
## skewhat standard error lower bound upper bound
## [1,] -0.3789 0.1669 -0.7127 -0.04505
kurt95CI(returns.df, "soybeans")
## kurthat standard error lower bound upper bound
## [1,] 2.014 0.4499 1.114 2.914
mu95CI(grains.z, "soymeal")
## muhat standard error lower bound upper bound
## [1,] 0.002199 0.001219 -0.0002378 0.004636
sd95CI(grains.z, "soymeal")
## sdhat standard error lower bound upper bound
## [1,] 0.03783 0.00117 0.03549 0.04017
skew95CI(returns.df, "soymeal")
## skewhat standard error lower bound upper bound
## [1,] -0.1422 0.1745 -0.4911 0.2068
kurt95CI(returns.df, "soymeal")
## kurthat standard error lower bound upper bound
## [1,] 1.81 0.5015 0.8072 2.813
head(rollMomentsGrains.ls[["soybeans"]][["mu"]])
## muhat.lower muhat muhat.upper
## 1993(41) 0.10080 0.15 0.1992
## 1993(42) 0.07043 0.12 0.1696
## 1993(43) 0.07043 0.12 0.1696
## 1993(44) 0.06064 0.11 0.1594
## 1993(45) 0.05068 0.10 0.1493
## 1993(46) 0.08008 0.13 0.1799
plotRollingMoments(data = rollMomentsGrains.ls, commo = "soybeans", statistic = "mu")
head(rollMomentsGrains.ls[["soybeans"]][["sd"]])
## sdhat.lower sdhat sdhat.upper
## 1993(41) 0.1703 0.1774 0.1845
## 1993(42) 0.1717 0.1787 0.1857
## 1993(43) 0.1719 0.1787 0.1855
## 1993(44) 0.1709 0.1780 0.1850
## 1993(45) 0.1706 0.1778 0.1851
## 1993(46) 0.1730 0.1800 0.1870
plotRollingMoments(data = rollMomentsGrains.ls, commo = "soybeans", statistic = "sd")
head(rollMomentsGrains.ls[["soybeans"]][["skew"]])
## skewhat.lower skewhat skewhat.upper
## 1993(41) -1.581 -0.3050 0.9706
## 1993(42) -1.497 -0.2439 1.0093
## 1993(43) -1.545 -0.2439 1.0568
## 1993(44) -1.474 -0.2300 1.0136
## 1993(45) -1.495 -0.2077 1.0791
## 1993(46) -1.459 -0.2421 0.9744
plotRollingMoments(data = rollMomentsGrains.ls, commo = "soybeans", statistic = "skew")
head(rollMomentsGrains.ls[["soybeans"]][["kurt"]])
## kurthat.lower kurthat kurthat.upper
## 1993(41) 0.02091 2.418 4.815
## 1993(42) 0.07039 2.250 4.429
## 1993(43) 0.16907 2.250 4.330
## 1993(44) 0.12356 2.324 4.524
## 1993(45) 0.14760 2.337 4.526
## 1993(46) 0.07688 2.134 4.191
plotRollingMoments(data = rollMomentsGrains.ls, commo = "soybeans", statistic = "kurt")
mu95CI(grains.z, "soyoil")
## muhat standard error lower bound upper bound
## [1,] 0.000332 0.00106 -0.001788 0.002452
sd95CI(grains.z, "soyoil")
## sdhat standard error lower bound upper bound
## [1,] 0.03291 0.0009685 0.03097 0.03485
skew95CI(returns.df, "soyoil")
## skewhat standard error lower bound upper bound
## [1,] -0.1591 0.1505 -0.4602 0.1419
kurt95CI(returns.df, "soyoil")
## kurthat standard error lower bound upper bound
## [1,] 1.397 0.3381 0.7209 2.073
head(rollMomentsGrains.ls[["soyoil"]][["mu"]])
## muhat.lower muhat muhat.upper
## 1993(41) 0.11735 0.17 0.2227
## 1993(42) 0.06677 0.12 0.1732
## 1993(43) 0.04701 0.10 0.1530
## 1993(44) 0.05681 0.11 0.1632
## 1993(45) 0.06663 0.12 0.1734
## 1993(46) 0.06663 0.12 0.1734
plotRollingMoments(data = rollMomentsGrains.ls, commo = "soyoil", statistic = "mu")
head(rollMomentsGrains.ls[["soyoil"]][["sd"]])
## sdhat.lower sdhat sdhat.upper
## 1993(41) 0.1846 0.1898 0.1950
## 1993(42) 0.1866 0.1919 0.1972
## 1993(43) 0.1859 0.1911 0.1962
## 1993(44) 0.1863 0.1918 0.1972
## 1993(45) 0.1873 0.1924 0.1976
## 1993(46) 0.1871 0.1924 0.1978
plotRollingMoments(data = rollMomentsGrains.ls, commo = "soyoil", statistic = "sd")
head(rollMomentsGrains.ls[["soyoil"]][["skew"]])
## skewhat.lower skewhat skewhat.upper
## 1993(41) -0.6462 -0.04066 0.5649
## 1993(42) -0.5770 0.02655 0.6301
## 1993(43) -0.5536 0.06497 0.6835
## 1993(44) -0.5284 0.04803 0.6245
## 1993(45) -0.5537 0.03114 0.6160
## 1993(46) -0.5499 0.03114 0.6122
plotRollingMoments(data = rollMomentsGrains.ls, commo = "soyoil", statistic = "skew")
head(rollMomentsGrains.ls[["soyoil"]][["kurt"]])
## kurthat.lower kurthat kurthat.upper
## 1993(41) -0.8054 0.14725 1.0999
## 1993(42) -0.9238 0.05458 1.0330
## 1993(43) -0.8813 0.10717 1.0956
## 1993(44) -0.8308 0.06498 0.9608
## 1993(45) -0.8653 0.02462 0.9146
## 1993(46) -0.8978 0.02462 0.9470
plotRollingMoments(data = rollMomentsGrains.ls, commo = "soyoil", statistic = "kurt")
see similar analysis for the other commodities classes here metals : http://rpubs.com/Bautheac/10353 energy : http://rpubs.com/Bautheac/10367 livestock : http://rpubs.com/Bautheac/10358 softs : http://rpubs.com/Bautheac/10355
see also preliminary workouts on correlations here http://rpubs.com/Bautheac/10352
See the functions called here https://gist.github.com/bautheac/b2e8269868a881b66a69/raw/