R version 2.15.1 (2012-06-22) – “Roasted Marshmallows”
Copyright © 2012 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit)
date()
## [1] "Tue Jul 10 20:57:50 2012"
Set working directory.
setwd("~/Dropbox/RicksLaw")
Get the required packages.
require(ggplot2)
require(quantreg)
require(R.matlab)
require(Renext)
require(grid)
Get additional functions. Power law fit.
source("plfit.R")
The data were sent via email from Rick Murnane on April 20, 2012 in a file called stormData.csv. Import data, remove missing damage and commas from losses and base amounts.
stormData = read.csv("stormData.csv", header = TRUE, na.strings = "NA")
stormData = stormData[!is.na(stormData$CURRENT.DAMAGE.2012), ]
stormData$Loss = as.numeric(gsub(",", "", stormData$CURRENT.DAMAGE.2012))
stormData$Base = as.numeric(gsub(",", "", stormData$BASE.DAMAGE))
Convert landfall wind speeds in mph to m/s and change date format.
stormData$W = stormData$WmaxLandfall * 0.44704
stormData$date = as.Date(stormData$format.change, format = "%b %d, %Y")
Plot a time series of the data.
p0 = ggplot(stormData, aes(x = date, y = Loss/1e+09)) + geom_point() +
theme_bw()
p0 = p0 + xlab("") + ylab("Economic Damage [billion USD (2012)]")
p0
First test to see whether a power law is an adequate description of these data.
plfit(stormData$Loss)
## $xmin
## [1] 3.54e+09
##
## $alpha
## [1] 1.856
##
## $D
## [1] 0.08856
##
A value of D < .1 indicates a power law is a poor fit for these data.
Next examine the evidence in support of an exponential fit. Start by generating random deviates from an exponential distribution. Set the number of deviates equal to the number of loss events.
n = length(stormData$Loss)
set.seed(1234)
x = gofExp.test(rexp(n))
x
## $statistic
## [1] 254
##
## $df
## [1] 237
##
## $p.value
## [1] 0.4275
##
## $method
## [1] "Bartlett gof for exponential"
##
As expected there is little evidence (large p-value) to reject the null hypothesis that the random deviates have an exponential distribution.
Now examine the loss data for losses exceeding $3.5 billion.
threshold = 3.5e+09
data = stormData$Loss[stormData$Loss > threshold]
x = gofExp.test(data)
x
## $statistic
## [1] 61.82
##
## $df
## [1] 66
##
## $p.value
## [1] 0.7539
##
## $method
## [1] "Bartlett gof for exponential"
##
Similar to the random data, there is no evidence to reject the null hypothesis that the loss data have an exponential distribution.
stormDataF = stormData[stormData$LANDFALL.STATE == "FL", ]
stormDataG = stormData[stormData$LANDFALL.STATE == "TX" | stormData$LANDFALL.STATE ==
"LA" | stormData$LANDFALL.STATE == "MS" | stormData$LANDFALL.STATE == "AL",
]
stormDataE = stormData[stormData$LANDFALL.STATE != "TX" & stormData$LANDFALL.STATE !=
"LA" & stormData$LANDFALL.STATE != "MS" & stormData$LANDFALL.STATE != "AL" &
stormData$LANDFALL.STATE != "FL", ]
Plot U.S. losses and losses by region on a semi-log plot.
p1 = ggplot(stormData, aes(x = W, y = Loss)) + geom_point() + scale_y_log10(limits = c(1e+06,
1e+12)) + theme_bw()
p1 = p1 + xlab("Wind speed [m/s]") + ylab("Economic Damage [billion USD (2012)]") +
xlim(10, 90)
p1 = p1 + stat_quantile(aes(colour = ..quantile..), quantiles = c(0.1,
0.25, 0.5, 0.75, 0.9))
p1 = p1 + scale_colour_gradient(low = "yellow", high = "red") + opts(legend.position = "none")
p2 = ggplot(stormDataG, aes(W, Loss)) + geom_point() + scale_y_log10(limits = c(1e+06,
1e+12)) + theme_bw()
p2 = p2 + xlab("Wind speed [m/s]") + ylab("Economic Damage [billion USD (2012)]") +
xlim(10, 90)
p2 = p2 + stat_quantile(aes(colour = ..quantile..), quantiles = c(0.1,
0.25, 0.5, 0.75, 0.9))
p2 = p2 + scale_colour_gradient(low = "yellow", high = "red") + opts(legend.position = "none")
p3 = ggplot(stormDataF, aes(W, Loss)) + geom_point() + scale_y_log10(limits = c(1e+06,
1e+12)) + theme_bw()
p3 = p3 + xlab("Wind speed [m/s]") + ylab("Economic Damage [billion USD (2012)]") +
xlim(10, 90)
p3 = p3 + stat_quantile(aes(colour = ..quantile..), quantiles = c(0.1,
0.25, 0.5, 0.75, 0.9))
p3 = p3 + scale_colour_gradient(low = "yellow", high = "red") + opts(legend.position = "none")
p4 = ggplot(stormDataE, aes(W, Loss)) + geom_point() + scale_y_log10(limits = c(1e+06,
1e+12)) + theme_bw()
p4 = p4 + xlab("Wind speed [m/s]") + ylab("Economic Damage [billion USD (2012)]") +
xlim(10, 90)
p4 = p4 + stat_quantile(aes(colour = ..quantile..), quantiles = c(0.1,
0.25, 0.5, 0.75, 0.9))
p4 = p4 + scale_colour_gradient(low = "yellow", high = "red") + opts(legend.position = "none")
p1 = p1 + opts(title = "U.S.", plot.title = theme_text(hjust = 0,
size = 15))
p2 = p2 + opts(title = "Gulf coast", plot.title = theme_text(hjust = 0,
size = 15))
p3 = p3 + opts(title = "Florida", plot.title = theme_text(hjust = 0,
size = 15))
p4 = p4 + opts(title = "East coast", plot.title = theme_text(hjust = 0,
size = 15))
grid.newpage()
pushViewport(viewport(layout = grid.layout(2, 2)))
vplayout = function(x, y) viewport(layout.pos.row = x, layout.pos.col = y)
print(p1, vp = vplayout(1, 1))
print(p2, vp = vplayout(1, 2))
print(p3, vp = vplayout(2, 1))
print(p4, vp = vplayout(2, 2))
summary(rq(log10(Loss) ~ I(W - mean(W)), data = stormData, tau = c(0.1,
0.25, 0.5, 0.75, 0.9)))
##
## Call: rq(formula = log10(Loss) ~ I(W - mean(W)), tau = c(0.1, 0.25,
## 0.5, 0.75, 0.9), data = stormData)
##
## tau: [1] 0.1
##
## Coefficients:
## coefficients lower bd upper bd
## (Intercept) 7.72155 7.46665 7.90240
## I(W - mean(W)) 0.05031 0.03867 0.06367
##
## Call: rq(formula = log10(Loss) ~ I(W - mean(W)), tau = c(0.1, 0.25,
## 0.5, 0.75, 0.9), data = stormData)
##
## tau: [1] 0.25
##
## Coefficients:
## coefficients lower bd upper bd
## (Intercept) 8.21825 8.15616 8.39654
## I(W - mean(W)) 0.05083 0.04504 0.05873
##
## Call: rq(formula = log10(Loss) ~ I(W - mean(W)), tau = c(0.1, 0.25,
## 0.5, 0.75, 0.9), data = stormData)
##
## tau: [1] 0.5
##
## Coefficients:
## coefficients lower bd upper bd
## (Intercept) 8.81986 8.73870 8.94752
## I(W - mean(W)) 0.04921 0.03666 0.05407
##
## Call: rq(formula = log10(Loss) ~ I(W - mean(W)), tau = c(0.1, 0.25,
## 0.5, 0.75, 0.9), data = stormData)
##
## tau: [1] 0.75
##
## Coefficients:
## coefficients lower bd upper bd
## (Intercept) 9.40036 9.29659 9.55117
## I(W - mean(W)) 0.04024 0.03222 0.05292
##
## Call: rq(formula = log10(Loss) ~ I(W - mean(W)), tau = c(0.1, 0.25,
## 0.5, 0.75, 0.9), data = stormData)
##
## tau: [1] 0.9
##
## Coefficients:
## coefficients lower bd upper bd
## (Intercept) 9.82012 9.74608 9.97771
## I(W - mean(W)) 0.03841 0.02772 0.04553
Interpretation. For every 1 m/s increase in landfall wind speed the industry losses increase by a factor of exp(.04921) or 1.050441 (5%).
Are the slopes significantly different across the quantiles? Use analysis of deviance.
fit0 = rq(log10(Loss) ~ I(W - mean(W)), data = stormData, tau = 0.1)
fit1 = rq(log10(Loss) ~ I(W - mean(W)), data = stormData, tau = 0.25)
fit2 = rq(log10(Loss) ~ I(W - mean(W)), data = stormData, tau = 0.5)
fit3 = rq(log10(Loss) ~ I(W - mean(W)), data = stormData, tau = 0.75)
fit4 = rq(log10(Loss) ~ I(W - mean(W)), data = stormData, tau = 0.9)
anova(fit0, fit1, fit2, fit3, fit4)
## Quantile Regression Analysis of Deviance Table
##
## Model: log10(Loss) ~ I(W - mean(W))
## Joint Test of Equality of Slopes: tau in { 0.1 0.25 0.5 0.75 0.9 }
##
## Df Resid Df F value Pr(>F)
## 1 4 1186 1.31 0.26
Interpretation. There is some evidence to reject the null hypothesis of equal slopes but not enough to reject it at a significance level of .1.
Repeat for winds greater than 20 m/s.
th = 20
fit0 = rq(log10(Loss) ~ I(W - mean(W)), data = stormData[stormData$W >
th, ], tau = 0.1)
fit1 = rq(log10(Loss) ~ I(W - mean(W)), data = stormData[stormData$W >
th, ], tau = 0.25)
fit2 = rq(log10(Loss) ~ I(W - mean(W)), data = stormData[stormData$W >
th, ], tau = 0.5)
fit3 = rq(log10(Loss) ~ I(W - mean(W)), data = stormData[stormData$W >
th, ], tau = 0.75)
fit4 = rq(log10(Loss) ~ I(W - mean(W)), data = stormData[stormData$W >
th, ], tau = 0.9)
anova(fit0, fit1, fit2, fit3, fit4)
## Quantile Regression Analysis of Deviance Table
##
## Model: log10(Loss) ~ I(W - mean(W))
## Joint Test of Equality of Slopes: tau in { 0.1 0.25 0.5 0.75 0.9 }
##
## Df Resid Df F value Pr(>F)
## 1 4 1126 2.82 0.024 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Repeat for winds greater than 35 m/s.
th = 35
fit0 = rq(log10(Loss) ~ I(W - mean(W)), data = stormData[stormData$W >
th, ], tau = 0.1)
fit1 = rq(log10(Loss) ~ I(W - mean(W)), data = stormData[stormData$W >
th, ], tau = 0.25)
fit2 = rq(log10(Loss) ~ I(W - mean(W)), data = stormData[stormData$W >
th, ], tau = 0.5)
fit3 = rq(log10(Loss) ~ I(W - mean(W)), data = stormData[stormData$W >
th, ], tau = 0.75)
fit4 = rq(log10(Loss) ~ I(W - mean(W)), data = stormData[stormData$W >
th, ], tau = 0.9)
anova(fit0, fit1, fit2, fit3, fit4)
## Warning: 1 non-positive fis
## Quantile Regression Analysis of Deviance Table
##
## Model: log10(Loss) ~ I(W - mean(W))
## Joint Test of Equality of Slopes: tau in { 0.1 0.25 0.5 0.75 0.9 }
##
## Df Resid Df F value Pr(>F)
## 1 4 751 1.54 0.19
Repeat for winds greater than 50 m/s.
th = 50
fit0 = rq(log10(Loss) ~ I(W - mean(W)), data = stormData[stormData$W >
th, ], tau = 0.1)
fit1 = rq(log10(Loss) ~ I(W - mean(W)), data = stormData[stormData$W >
th, ], tau = 0.25)
fit2 = rq(log10(Loss) ~ I(W - mean(W)), data = stormData[stormData$W >
th, ], tau = 0.5)
## Warning: Solution may be nonunique
fit3 = rq(log10(Loss) ~ I(W - mean(W)), data = stormData[stormData$W >
th, ], tau = 0.75)
fit4 = rq(log10(Loss) ~ I(W - mean(W)), data = stormData[stormData$W >
th, ], tau = 0.9)
anova(fit0, fit1, fit2, fit3, fit4)
## Warning: 1 non-positive fis
## Quantile Regression Analysis of Deviance Table
##
## Model: log10(Loss) ~ I(W - mean(W))
## Joint Test of Equality of Slopes: tau in { 0.1 0.25 0.5 0.75 0.9 }
##
## Df Resid Df F value Pr(>F)
## 1 4 301 0.59 0.67
Repeat for Gulf coast losses.
summary(rq(log10(Loss) ~ I(W - mean(W)), data = stormDataG, tau = c(0.1,
0.25, 0.5, 0.75, 0.9)))
##
## Call: rq(formula = log10(Loss) ~ I(W - mean(W)), tau = c(0.1, 0.25,
## 0.5, 0.75, 0.9), data = stormDataG)
##
## tau: [1] 0.1
##
## Coefficients:
## coefficients lower bd upper bd
## (Intercept) 7.74467 7.40321 8.02270
## I(W - mean(W)) 0.05031 0.03578 0.06409
##
## Call: rq(formula = log10(Loss) ~ I(W - mean(W)), tau = c(0.1, 0.25,
## 0.5, 0.75, 0.9), data = stormDataG)
##
## tau: [1] 0.25
##
## Coefficients:
## coefficients lower bd upper bd
## (Intercept) 8.23030 8.10343 8.39001
## I(W - mean(W)) 0.05118 0.04688 0.05889
##
## Call: rq(formula = log10(Loss) ~ I(W - mean(W)), tau = c(0.1, 0.25,
## 0.5, 0.75, 0.9), data = stormDataG)
##
## tau: [1] 0.5
##
## Coefficients:
## coefficients lower bd upper bd
## (Intercept) 8.82516 8.63912 8.94991
## I(W - mean(W)) 0.04414 0.03748 0.05699
##
## Call: rq(formula = log10(Loss) ~ I(W - mean(W)), tau = c(0.1, 0.25,
## 0.5, 0.75, 0.9), data = stormDataG)
##
## tau: [1] 0.75
##
## Coefficients:
## coefficients lower bd upper bd
## (Intercept) 9.24373 9.12324 9.59427
## I(W - mean(W)) 0.05005 0.02957 0.06549
##
## Call: rq(formula = log10(Loss) ~ I(W - mean(W)), tau = c(0.1, 0.25,
## 0.5, 0.75, 0.9), data = stormDataG)
##
## tau: [1] 0.9
##
## Coefficients:
## coefficients lower bd upper bd
## (Intercept) 9.75462 9.67345 9.99286
## I(W - mean(W)) 0.03974 0.02659 0.05883
fit0 = rq(log10(Loss) ~ I(W - mean(W)), data = stormDataG, tau = 0.1)
fit1 = rq(log10(Loss) ~ I(W - mean(W)), data = stormDataG, tau = 0.25)
fit2 = rq(log10(Loss) ~ I(W - mean(W)), data = stormDataG, tau = 0.5)
fit3 = rq(log10(Loss) ~ I(W - mean(W)), data = stormDataG, tau = 0.75)
fit4 = rq(log10(Loss) ~ I(W - mean(W)), data = stormDataG, tau = 0.9)
anova(fit0, fit1, fit2, fit3, fit4)
## Quantile Regression Analysis of Deviance Table
##
## Model: log10(Loss) ~ I(W - mean(W))
## Joint Test of Equality of Slopes: tau in { 0.1 0.25 0.5 0.75 0.9 }
##
## Df Resid Df F value Pr(>F)
## 1 4 481 0.87 0.48
Repeat for Florida losses.
summary(rq(log10(Loss) ~ I(W - mean(W)), data = stormDataF, tau = c(0.1,
0.25, 0.5, 0.75, 0.9)))
##
## Call: rq(formula = log10(Loss) ~ I(W - mean(W)), tau = c(0.1, 0.25,
## 0.5, 0.75, 0.9), data = stormDataF)
##
## tau: [1] 0.1
##
## Coefficients:
## coefficients lower bd upper bd
## (Intercept) 7.85737 7.46145 7.96795
## I(W - mean(W)) 0.04807 0.03494 0.06806
##
## Call: rq(formula = log10(Loss) ~ I(W - mean(W)), tau = c(0.1, 0.25,
## 0.5, 0.75, 0.9), data = stormDataF)
##
## tau: [1] 0.25
##
## Coefficients:
## coefficients lower bd upper bd
## (Intercept) 8.30576 8.10463 8.60113
## I(W - mean(W)) 0.05539 0.03295 0.06574
##
## Call: rq(formula = log10(Loss) ~ I(W - mean(W)), tau = c(0.1, 0.25,
## 0.5, 0.75, 0.9), data = stormDataF)
##
## tau: [1] 0.5
##
## Coefficients:
## coefficients lower bd upper bd
## (Intercept) 9.01224 8.79746 9.08796
## I(W - mean(W)) 0.04510 0.03738 0.05220
##
## Call: rq(formula = log10(Loss) ~ I(W - mean(W)), tau = c(0.1, 0.25,
## 0.5, 0.75, 0.9), data = stormDataF)
##
## tau: [1] 0.75
##
## Coefficients:
## coefficients lower bd upper bd
## (Intercept) 9.41112 9.29407 9.65128
## I(W - mean(W)) 0.03955 0.03430 0.05382
##
## Call: rq(formula = log10(Loss) ~ I(W - mean(W)), tau = c(0.1, 0.25,
## 0.5, 0.75, 0.9), data = stormDataF)
##
## tau: [1] 0.9
##
## Coefficients:
## coefficients lower bd upper bd
## (Intercept) 9.83363 9.69133 10.10720
## I(W - mean(W)) 0.04069 0.02485 0.04971
fit0 = rq(log10(Loss) ~ I(W - mean(W)), data = stormDataF, tau = 0.1)
fit1 = rq(log10(Loss) ~ I(W - mean(W)), data = stormDataF, tau = 0.25)
fit2 = rq(log10(Loss) ~ I(W - mean(W)), data = stormDataF, tau = 0.5)
fit3 = rq(log10(Loss) ~ I(W - mean(W)), data = stormDataF, tau = 0.75)
fit4 = rq(log10(Loss) ~ I(W - mean(W)), data = stormDataF, tau = 0.9)
anova(fit0, fit1, fit2, fit3, fit4)
## Quantile Regression Analysis of Deviance Table
##
## Model: log10(Loss) ~ I(W - mean(W))
## Joint Test of Equality of Slopes: tau in { 0.1 0.25 0.5 0.75 0.9 }
##
## Df Resid Df F value Pr(>F)
## 1 4 406 0.59 0.67
Repeat for East coast losses.
summary(rq(log10(Loss) ~ I(W - mean(W)), data = stormDataE, tau = c(0.1,
0.25, 0.5, 0.75, 0.9)))
##
## Call: rq(formula = log10(Loss) ~ I(W - mean(W)), tau = c(0.1, 0.25,
## 0.5, 0.75, 0.9), data = stormDataE)
##
## tau: [1] 0.1
##
## Coefficients:
## coefficients lower bd upper bd
## (Intercept) 7.33160 6.68642 7.93765
## I(W - mean(W)) 0.05655 0.02163 0.11378
##
## Call: rq(formula = log10(Loss) ~ I(W - mean(W)), tau = c(0.1, 0.25,
## 0.5, 0.75, 0.9), data = stormDataE)
##
## tau: [1] 0.25
##
## Coefficients:
## coefficients lower bd upper bd
## (Intercept) 8.18780 7.96350 8.42834
## I(W - mean(W)) 0.04762 0.02888 0.07270
##
## Call: rq(formula = log10(Loss) ~ I(W - mean(W)), tau = c(0.1, 0.25,
## 0.5, 0.75, 0.9), data = stormDataE)
##
## tau: [1] 0.5
##
## Coefficients:
## coefficients lower bd upper bd
## (Intercept) 8.72353 8.51642 8.94263
## I(W - mean(W)) 0.05022 0.01212 0.08331
##
## Call: rq(formula = log10(Loss) ~ I(W - mean(W)), tau = c(0.1, 0.25,
## 0.5, 0.75, 0.9), data = stormDataE)
##
## tau: [1] 0.75
##
## Coefficients:
## coefficients lower bd upper bd
## (Intercept) 9.46691 9.33290 9.76235
## I(W - mean(W)) 0.03462 0.00653 0.05574
##
## Call: rq(formula = log10(Loss) ~ I(W - mean(W)), tau = c(0.1, 0.25,
## 0.5, 0.75, 0.9), data = stormDataE)
##
## tau: [1] 0.9
##
## Coefficients:
## coefficients lower bd upper bd
## (Intercept) 10.05263 9.76661 10.31862
## I(W - mean(W)) 0.01537 -0.01226 0.06381
fit0 = rq(log10(Loss) ~ I(W - mean(W)), data = stormDataE, tau = 0.1)
fit1 = rq(log10(Loss) ~ I(W - mean(W)), data = stormDataE, tau = 0.25)
fit2 = rq(log10(Loss) ~ I(W - mean(W)), data = stormDataE, tau = 0.5)
fit3 = rq(log10(Loss) ~ I(W - mean(W)), data = stormDataE, tau = 0.75)
fit4 = rq(log10(Loss) ~ I(W - mean(W)), data = stormDataE, tau = 0.9)
anova(fit0, fit1, fit2, fit3, fit4)
## Quantile Regression Analysis of Deviance Table
##
## Model: log10(Loss) ~ I(W - mean(W))
## Joint Test of Equality of Slopes: tau in { 0.1 0.25 0.5 0.75 0.9 }
##
## Df Resid Df F value Pr(>F)
## 1 4 291 0.66 0.62
Plot of the quantile regression losses on landfall wind speed.
modelUS = rq(log10(Loss) ~ I(W - mean(W)), data = stormData, tau = seq(0.1,
0.9, 0.05))
modelG = rq(log10(Loss) ~ I(W - mean(W)), data = stormDataG, tau = seq(0.1,
0.9, 0.05))
modelF = rq(log10(Loss) ~ I(W - mean(W)), data = stormDataF, tau = seq(0.1,
0.9, 0.05))
modelE = rq(log10(Loss) ~ I(W - mean(W)), data = stormDataE, tau = seq(0.1,
0.9, 0.05))
par(mfrow = c(2, 2))
plot(summary(modelUS, se = "boot"), parm = 2, xaxt = "n", mar = c(5,
5, 4, 2) + 0.1, pch = 16, lwd = 2, ylab = "Slope", xlab = "Wind speed quantile",
main = "U.S. Losses")
plot(summary(modelG, se = "boot"), parm = 2, xaxt = "n", mar = c(5,
5, 4, 2) + 0.1, pch = 16, lwd = 2, ylab = "Slope", xlab = "Wind speed quantile",
main = "Gulf coast Losses")
plot(summary(modelF, se = "boot"), parm = 2, xaxt = "n", mar = c(5,
5, 4, 2) + 0.1, pch = 16, lwd = 2, ylab = "Slope", xlab = "Wind speed quantile",
main = "Florida Losses")
plot(summary(modelE, se = "boot"), parm = 2, xaxt = "n", mar = c(5,
5, 4, 2) + 0.1, pch = 16, lwd = 2, ylab = "Slope", xlab = "Wind speed quantile",
main = "East coast Losses")
There is no statistically significant change in the slope over the quantiles of wind speed. This result is consistent with the anova test above.
Repeat using a power-law model.
summary(rq(log10(Loss) ~ log10(W), data = stormData, tau = c(0.1,
0.25, 0.5, 0.75, 0.9)))
##
## Call: rq(formula = log10(Loss) ~ log10(W), tau = c(0.1, 0.25, 0.5,
## 0.75, 0.9), data = stormData)
##
## tau: [1] 0.1
##
## Coefficients:
## coefficients lower bd upper bd
## (Intercept) -0.2497 -1.4929 3.0785
## log10(W) 4.9828 2.6884 5.8636
##
## Call: rq(formula = log10(Loss) ~ log10(W), tau = c(0.1, 0.25, 0.5,
## 0.75, 0.9), data = stormData)
##
## tau: [1] 0.25
##
## Coefficients:
## coefficients lower bd upper bd
## (Intercept) 0.7744 -0.7908 2.2920
## log10(W) 4.7511 4.0017 5.6343
##
## Call: rq(formula = log10(Loss) ~ log10(W), tau = c(0.1, 0.25, 0.5,
## 0.75, 0.9), data = stormData)
##
## tau: [1] 0.5
##
## Coefficients:
## coefficients lower bd upper bd
## (Intercept) 2.6810 0.8291 3.6395
## log10(W) 3.8989 3.2829 4.9993
##
## Call: rq(formula = log10(Loss) ~ log10(W), tau = c(0.1, 0.25, 0.5,
## 0.75, 0.9), data = stormData)
##
## tau: [1] 0.75
##
## Coefficients:
## coefficients lower bd upper bd
## (Intercept) 3.749 2.937 4.725
## log10(W) 3.579 3.010 4.056
##
## Call: rq(formula = log10(Loss) ~ log10(W), tau = c(0.1, 0.25, 0.5,
## 0.75, 0.9), data = stormData)
##
## tau: [1] 0.9
##
## Coefficients:
## coefficients lower bd upper bd
## (Intercept) 4.901 4.030 6.086
## log10(W) 3.105 2.358 3.744
fit0 = rq(log10(Loss) ~ log10(W), data = stormData, tau = 0.1)
fit1 = rq(log10(Loss) ~ log10(W), data = stormData, tau = 0.25)
fit2 = rq(log10(Loss) ~ log10(W), data = stormData, tau = 0.5)
fit3 = rq(log10(Loss) ~ log10(W), data = stormData, tau = 0.75)
fit4 = rq(log10(Loss) ~ log10(W), data = stormData, tau = 0.9)
anova(fit0, fit1, fit2, fit3, fit4)
## Quantile Regression Analysis of Deviance Table
##
## Model: log10(Loss) ~ log10(W)
## Joint Test of Equality of Slopes: tau in { 0.1 0.25 0.5 0.75 0.9 }
##
## Df Resid Df F value Pr(>F)
## 1 4 1186 2.01 0.091 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Plot of the quantile regression losses on landfall wind speed.
modelUS = rq(log10(Loss) ~ log10(W), data = stormData, tau = seq(0.1,
0.9, 0.05))
plot(summary(modelUS, se = "boot"), parm = 2, xaxt = "n", mar = c(5,
5, 4, 2) + 0.1, pch = 16, lwd = 2, ylab = "Scaling exponent", xlab = "Wind speed quantile",
main = "U.S. Losses (power law)")
fit0I = rq(log10(Loss) ~ W, data = stormDataE, tau = 0.1)
fit2I = rq(log10(Loss) ~ W, data = stormDataE, tau = 0.5)
fit4I = rq(log10(Loss) ~ W, data = stormDataE, tau = 0.9)
fit5I = rq(log10(Loss) ~ W, data = stormDataE, tau = 0.95)
10^predict(fit4I, newdata = data.frame(W = c(39, 28)), interval = "confidence",
level = 0.9)
## fit lower higher
## 1 1.106e+10 5.110e+09 2.392e+10
## 2 7.489e+09 2.035e+09 2.756e+10