<p><div><div><div><div data-width="100" data-show-faces="false" data-send="true" data-layout="button_count" data-height="21" data-href="https://www.r-bloggers.com/forcasting-natural-catastrophes-is-rather-difficult/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://freakonometrics.hypotheses.org/50057"> R-english – Freakonometrics</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>Following my <a target="_blank" rel="nofollow" href="https://freakonometrics.hypotheses.org/50003">previous post</a>, I wanted to spend more time, on the time series with “global weather-related disaster losses as a proportion of global GDP” over the time period 1990-2016 that Roger Pilke sent me last night.</p>
<p><code>db=data.frame(year=1990:2016,<br>
ratio=c(.23,.27,.32,.37,.22,.26,.29,.15,.40,.28,.14,.09,.24,.18,.29,.51,.13,.17,.25,.13,.21,.29,.25,.2,.15,.12,.12))</code></p>
<p>In my <a target="_blank" rel="nofollow" href="https://freakonometrics.hypotheses.org/50003">previous post</a>, I spend some time explaining that we should provide some sort of ‘confidence interval’ when we try to predict a pattern. That was what we call ‘model uncertainty’. But there are two (important) issues that I did not mention. (1) it is a time series, so why not use techniques dedicated to time series objects ? (2) we do not really care actually about ‘model uncertainty’ (unless we want to assess if a decreasing trend is significant, or not), and we care more about real prediction uncertainty: <em>in the next ten years, what could be the range for the this ratio, with some given probability (say 95%)</em>? Could we say that with 95% chance the global weather-related disaster losses as a proportion of global GDP should be (each year) within 0 and 0.35 or 0 and 0.7?</p>
<p>A first idea might be to use exponential smoothing techniques (without a seasonal component here).</p>
<p><code>ratio=ts(db$ratio,start=1990,frequency=1)<br>
plot(ratio,xlim=c(1990,2030))<br>
hw=HoltWinters(ratio,gamma=FALSE)<br>
phw=predict(hw,n.ahead=15,prediction.interval = TRUE)<br>
plot(hw,phw,xlim=c(1990,2030))<br>
polygon(c(2017:2031,rev(2017:2031)), c(phw[,2],rev(phw[,3])),border=NA,col=rgb(0,0,1,.2))</code></p>
<p><img alt="" src="https://i1.wp.com/freakonometrics.hypotheses.org/files/2017/01/pilke11.png?w=456&ssl=1" data-recalc-dims="1"></p>
<p>The decreasing trend is coming from the fact that exponential smoothing is here a linear regression, with weight exponentially decaying with time (the older, the smaller the weight). But we cannot use that prediction, since the ratio cannot (obviously) be negative. So why not consider, here, the logarithm of the ratio</p>
<p><code>plot(db$year,db$ratio,type=&quot;l&quot;,xlim=c(1990,2030),ylim=c(-.2,.7),xlab=&quot;year&quot;,ylab=&quot;ratio&quot;)<br>
hw=HoltWinters(log(ratio),gamma=FALSE)<br>
phw=predict(hw,n.ahead=15,prediction.interval = TRUE)<br>
abline(v=2016,lty=2,col=&quot;grey&quot;)<br>
lines(2017:2031,exp(phw[,2]),col=&quot;blue&quot;)<br>
lines(2017:2031,exp(phw[,3]),col=&quot;blue&quot;)<br>
lines(c(1992:2016,2017:2031),c(exp(hw$fitted[,1]),exp(phw[,1])),col=&quot;red&quot;)<br>
polygon(c(2017:2031,rev(2017:2031)),exp(c(phw[,2],rev(phw[,3]))),border=NA,col=rgb(0,0,1,.2))<br>
abline(h=0,lty=2)</code></p>
<p><img alt="" src="https://i2.wp.com/freakonometrics.hypotheses.org/files/2017/01/pilke12.png?w=456&ssl=1" data-recalc-dims="1"></p>
<p>The confidence band is huge, here. What if we consider some ARIMA model here?</p>
<p><code>fit=auto.arima(ratio)<br>
farma=forecast(fit,15)<br>
farma=cbind(as.numeric(farma$fitted)[1:15],as.numeric(farma$lower[,1]),as.numeric(farma$upper[,1]),as.numeric(farma$lower[,2]),as.numeric(farma$upper[,2]))<br>
plot(db$year,db$ratio,type=&quot;l&quot;,xlim=c(1990,2030),ylim=c(-.2,.7),xlab=&quot;year&quot;,ylab=&quot;ratio&quot;)<br>
abline(v=2016,lty=2,col=&quot;grey&quot;)<br>
lines(2017:2031,farma[,4],col=&quot;blue&quot;)<br>
lines(2017:2031,farma[,5],col=&quot;blue&quot;)<br>
lines(2017:2031,farma[,1],col=&quot;red&quot;)<br>
polygon(c(2017:2031,rev(2017:2031)),c(farma[,4],rev(farma[,5])),border=NA,col=rgb(0,0,1,.2))<br>
abline(h=0,lty=2)</code></p>
<p><img alt="" src="https://i2.wp.com/freakonometrics.hypotheses.org/files/2017/01/pilke13.png?w=456&ssl=1" data-recalc-dims="1"></p>
<p>Here, there is an intercept, but no dynamics for the time series (which is considered, here, as a pure white noise). We get exactly the same if we consider the average value of the series</p>
<p><code>fit=lm(ratio~1,data=db)<br>
s=summary(fit)$sigma<br>
plot(db$year,db$ratio,type=&quot;l&quot;,xlim=c(1990,2030),ylim=c(-.2,.7),xlab=&quot;year&quot;,ylab=&quot;ratio&quot;)<br>
abline(v=2016,lty=2,col=&quot;grey&quot;)<br>
ndb=data.frame(year=2017:2031)<br>
pf=predict(fit,newdata=ndb)<br>
farma=cbind(pf,pf-1.96*s,pf+1.96*s)<br>
lines(2017:2031,farma[,2],col=&quot;blue&quot;)<br>
lines(2017:2031,farma[,3],col=&quot;blue&quot;)<br>
lines(1990:2031,c(predict(fit),farma[,1]),col=&quot;red&quot;)<br>
polygon(c(2017:2031,rev(2017:2031)),c(farma[,2],rev(farma[,3])),border=NA,col=rgb(0,0,1,.2))<br>
abline(h=0,lty=2)</code></p>
<p><img alt="" src="https://i2.wp.com/freakonometrics.hypotheses.org/files/2017/01/pilke19.png?w=456&ssl=1" data-recalc-dims="1"></p>
<p>Here, we get back to my <a target="_blank" rel="nofollow" href="https://freakonometrics.hypotheses.org/50003">previous post</a>, if we want to consider a possible trend (and not only an intercept)</p>
<p><code>fit=lm(ratio~year,data=db)<br>
s=summary(fit)$sigma<br>
plot(db$year,db$ratio,type=&quot;l&quot;,xlim=c(1990,2030),ylim=c(-.2,.7),xlab=&quot;year&quot;,ylab=&quot;ratio&quot;)<br>
abline(v=2016,lty=2,col=&quot;grey&quot;)<br>
ndb=data.frame(year=2017:2031)<br>
pf=predict(fit,newdata=ndb)<br>
farma=cbind(pf,pf-1.96*s,pf+1.96*s)<br>
lines(2017:2031,farma[,2],col=&quot;blue&quot;)<br>
lines(2017:2031,farma[,3],col=&quot;blue&quot;)<br>
lines(1990:2031,c(predict(fit),farma[,1]),col=&quot;red&quot;)<br>
polygon(c(2017:2031,rev(2017:2031)),c(farma[,2],rev(farma[,3])),border=NA,col=rgb(0,0,1,.2))<br>
abline(h=0,lty=2)</code></p>
<p><img alt="" src="https://i1.wp.com/freakonometrics.hypotheses.org/files/2017/01/pilke16.png?w=456&ssl=1" data-recalc-dims="1"></p>
<p>Again, the confidence region is not based on inference related error, but on model uncertainty: we try to visualize <em>where</em> future observations might be with (say) 95% chance. Note we can also consider (why not?) a quadratic regression</p>
<p><code>fit=lm(ratio~poly(year,2),data=db)<br>
s=summary(fit)$sigma<br>
plot(db$year,db$ratio,type=&quot;l&quot;,xlim=c(1990,2030),ylim=c(-.2,.7),xlab=&quot;year&quot;,ylab=&quot;ratio&quot;)<br>
abline(v=2016,lty=2,col=&quot;grey&quot;)<br>
ndb=data.frame(year=2017:2031)<br>
pf=predict(fit,newdata=ndb)<br>
farma=cbind(pf,pf-1.96*s,pf+1.96*s)<br>
lines(2017:2031,farma[,2],col=&quot;blue&quot;)<br>
lines(2017:2031,farma[,3],col=&quot;blue&quot;)<br>
lines(1990:2031,c(predict(fit),farma[,1]),col=&quot;red&quot;)<br>
polygon(c(2017:2031,rev(2017:2031)),c(farma[,2],rev(farma[,3])),border=NA,col=rgb(0,0,1,.2))<br>
abline(h=0,lty=2)</code></p>
<p><img alt="" src="https://i0.wp.com/freakonometrics.hypotheses.org/files/2017/01/pilke17.png?w=456&ssl=1" data-recalc-dims="1"></p>
<p>I am usually not a huge fan of those polynomial regression, but recently, I’ve seen that a lot in economic papers (like “if it’s not linear, add a squared version of the explanatory variable”, which is a rather odd strategy, I’ll publish some posts on that issue this year).</p>
<p>Here again, it might be more clever to consider a logarithmic transformation of the ratio, to insure that the ratio remains positive</p>
<p><code>fit=lm(log(ratio)~year,data=db)<br>
s=summary(fit)$sigma<br>
plot(db$year,db$ratio,type=&quot;l&quot;,xlim=c(1990,2030),ylim=c(-.2,.7),xlab=&quot;year&quot;,ylab=&quot;ratio&quot;)<br>
abline(v=2016,lty=2,col=&quot;grey&quot;)<br>
ndb=data.frame(year=2017:2031)<br>
pf=predict(fit,newdata=ndb)<br>
farma=cbind(exp(pf+s^2/2),exp(pf-1.96*s),exp(pf+1.96*s))<br>
lines(2017:2031,farma[,2],col=&quot;blue&quot;)<br>
lines(2017:2031,farma[,3],col=&quot;blue&quot;)<br>
lines(1990:2031,c(exp(predict(fit)+s^2/2),farma[,1]),col=&quot;red&quot;)<br>
polygon(c(2017:2031,rev(2017:2031)),c(farma[,2],rev(farma[,3])),border=NA,col=rgb(0,0,1,.2))<br>
abline(h=0,lty=2)</code></p>
<p><img alt="" src="https://i0.wp.com/freakonometrics.hypotheses.org/files/2017/01/pilke18.png?w=456&ssl=1" data-recalc-dims="1"></p>
<p>Observe that future trend is mainly driven by the three latest observations, that were rather low (compared with older observations). What if we remove them?</p>
<p><code>dbna=db<br>
db$ratio[25:27]=NA<br>
fit=lm(ratio~1,data=dbna)<br>
s=summary(fit)$sigma<br>
plot(db$year,db$ratio,type=&quot;l&quot;,xlim=c(1990,2030),ylim=c(-.2,.7),xlab=&quot;year&quot;,ylab=&quot;ratio&quot;)<br>
abline(v=2016-3,lty=2,col=&quot;grey&quot;)<br>
ndb=data.frame(year=2014:2031)<br>
pf=predict(fit,newdata=ndb)<br>
farma=cbind(pf,pf-1.96*s,pf+1.96*s)<br>
lines(2014:2031,farma[,2],col=&quot;blue&quot;)<br>
lines(2014:2031,farma[,3],col=&quot;blue&quot;)<br>
lines(1990:2031,c(predict(fit)[1:24],farma[,1]),col=&quot;red&quot;)<br>
polygon(c(2014:2031,rev(2014:2031)),c(farma[,2],rev(farma[,3])),border=NA,col=rgb(0,0,1,.2))<br>
abline(h=0,lty=2)</code></p>
<p><img alt="" src="https://i2.wp.com/freakonometrics.hypotheses.org/files/2017/01/pilke20.png?w=456&ssl=1" data-recalc-dims="1"></p>
<p>More funny, if we consider a quadratic regression, we obtain an increasing trend for the future</p>
<p><code>fit=lm(ratio~poly(year,2),data=dbna)<br>
s=summary(fit)$sigma<br>
plot(db$year,db$ratio,type=&quot;l&quot;,xlim=c(1990,2030),ylim=c(-.2,.7),xlab=&quot;year&quot;,ylab=&quot;ratio&quot;)<br>
abline(v=2016-3,lty=2,col=&quot;grey&quot;)<br>
ndb=data.frame(year=2014:2031)<br>
pf=predict(fit,newdata=ndb)<br>
farma=cbind(pf,pf-1.96*s,pf+1.96*s)<br>
lines(2014:2031,farma[,2],col=&quot;blue&quot;)<br>
lines(2014:2031,farma[,3],col=&quot;blue&quot;)<br>
lines(1990:2031,c(predict(fit)[1:24],farma[,1]),col=&quot;red&quot;)<br>
polygon(c(2014:2031,rev(2014:2031)),c(farma[,2],rev(farma[,3])),border=NA,col=rgb(0,0,1,.2))<br>
abline(h=0,lty=2)</code></p>
<p><img alt="" src="https://i0.wp.com/freakonometrics.hypotheses.org/files/2017/01/pilke22.png?w=456&ssl=1" data-recalc-dims="1"></p>
<p>As we can see, it is rather difficult to get relevant prediction for the future, based on 25 observations…. If anyone has a suggestion, comments are open…</p>
<p> </p>
<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="http://freakonometrics.hypotheses.org/50057"> R-english – Freakonometrics</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a title="The R Project for Statistical Computing" rel="nofollow" href="https://www.r-project.org/">R</a> news and <a title="R tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/tutorial">tutorials</a> on topics such as: <a title="Data science" rel="nofollow" href="https://www.r-bloggers.com/search/data%20science">Data science</a>, <a title="Big Data" rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data">Big Data, <a title="R jobs" rel="nofollow" href="https://www.r-users.com/">R jobs</a>, visualization (<a title="ggplot and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2">ggplot2</a>, <a title="Boxplots using lattice and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/boxplot">Boxplots</a>, <a title="Maps and gis" rel="nofollow" href="https://www.r-bloggers.com/search/map">maps</a>, <a title="Animation in R" rel="nofollow" href="https://www.r-bloggers.com/search/animation">animation</a>), programming (<a title="RStudio IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/RStudio">RStudio</a>, <a title="Sweave and literate programming" rel="nofollow" href="https://www.r-bloggers.com/search/sweave">Sweave</a>, <a title="LaTeX in R" rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX">LaTeX</a>, <a title="SQL and databases" rel="nofollow" href="https://www.r-bloggers.com/search/SQL">SQL</a>, <a title="Eclipse IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/eclipse">Eclipse</a>, <a title="git and github, Version Control System" rel="nofollow" href="https://www.r-bloggers.com/search/git">git</a>, <a title="Large data in R using Hadoop" rel="nofollow" href="https://www.r-bloggers.com/search/hadoop">hadoop</a>, <a title="Web Scraping of google, facebook, yahoo, twitter and more using R" rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping">Web Scraping</a>) statistics (<a title="Regressions and ANOVA analysis tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/regression">regression</a>, <a title="principal component analysis tutorial" rel="nofollow" href="https://www.r-bloggers.com/search/PCA">PCA</a>, <a title="Time series" rel="nofollow" href="https://www.r-bloggers.com/search/time+series">time series</a>, <a title="finance trading" rel="nofollow" href="https://www.r-bloggers.com/search/trading">trading</a>) and more...
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/6osYhR_FBPY"></a></div>