机器学习实践秘籍-四-

机器学习实践秘籍（四）

原文：annas-archive.org/md5/15b58dbdde87ea8cb4894e5d29bf0db5

译者：飞龙

协议：CC BY-NC-SA 4.0

第十章：神经网络

在本章中，我们将介绍以下内容：

• 模拟标准普尔 500

• 测量失业率

简介

神经网络：神经网络是一个有序的三元组 ，其中 是神经元的集合， 是一个集合 ，其元素被称为神经元 和神经元 之间的连接。函数 定义了权重，其中 是神经元 和神经元 之间连接的权重。数据通过连接在神经元之间传输，连接权重可以是兴奋的或抑制的。

模拟标准普尔 500

根据市值，纽约证券交易所或纳斯达克综合指数上市的 500 家最大公司的股票市值是通过标准普尔 500 指数来衡量的。标准普尔提供了一种基于股票价格的市场和经济走势的快速观察。标准普尔 500 指数是金融媒体和专业人士最常用的衡量指标。标准普尔 500 指数是通过将所有标准普尔 500 股票的调整市值相加，然后除以标准普尔开发的指数除数来计算的。当出现股票分割、特殊股息或分拆等可能影响指数价值的情况时，除数会进行调整。除数确保这些非经济因素不会影响指数。

准备工作

为了使用神经网络来模拟标准普尔 500 指数，我们将使用从 GSPC 数据集收集的数据集。

第 1 步 - 收集和描述数据

要使用的数据集是 2009 年 1 月 1 日至 2014 年 1 月 1 日之间的 GSPC 每日收盘股票价值。这个数据集在 www.yahoo.com/ 上免费提供，我们将从那里下载数据。

如何操作...

让我们深入了解细节。

第 2 步 - 探索数据

首先，需要加载以下包：

 > install.packages("quantmod")
 > install.packages("neuralnet")
 > library(quantmod)
 > library(neuralnet)

让我们下载数据。我们将首先标记所需时间段的开始和结束日期。

as.Date() 函数用于将字符表示和 Date 类对象转换为日历日期。

数据集的开始日期存储在 startDate 中，它代表日历日期的字符向量表示。这种表示的格式是 YYYY-MM-DD：

> startDate = as.Date("2009-01-01")

数据集的结束日期存储在 endDate 中，它代表日历日期的字符向量表示。这种表示的格式是 YYYY-MM-DD：

> endDate = as.Date("2014-01-01")

使用getSymbols()函数加载数据：该函数从多个来源加载数据，无论是本地还是远程。GSPC是字符向量，指定要加载的符号名称。src = yahoo指定了数据来源方法：

> getSymbols("^GSPC", src="img/yahoo", from=startDate, to=endDate)

步骤 3 - 计算指标

计算相对强弱指数：这是最近上升价格变动与绝对价格变动的比率。使用RSI()函数计算相对强弱指数。GSPC数据框用作价格序列。n = 3代表移动平均的周期数。结果存储在relativeStrengthIndex3数据框中：

> relativeStrengthIndex3 <- RSI(Op(GSPC),n=3)

探索价格变化的总结：为此使用summary()函数。该函数提供一系列描述性统计，以生成relativeStrengthIndex3数据框的结果摘要：

> summary(relativeStrengthIndex3)

结果如下：

EMA()函数使用GSPC符号作为价格序列。n = 5代表平均的时间周期。结果存储在exponentialMovingAverage5数据框中：

> exponentialMovingAverage5 <- EMA(Op(GSPC),n=5)

打印exponentialMovingAverage5数据框：head()函数返回exponentialMovingAverage5数据框的前部分。exponentialMovingAverage5数据框作为输入参数传递：

> head(exponentialMovingAverage5)

结果如下：

探索价格变化的总结。为此，使用summary()函数。此函数提供一系列描述性统计，以生成exponentialMovingAverage5数据框的结果摘要。

> summary(exponentialMovingAverage5)

结果如下：

计算GSPC和exponentialMovingAverage5的指数开盘价之间的差异：

> exponentialMovingAverageDiff <- Op(GSPC) - exponentialMovingAverage5

现在让我们打印exponentialMovingAverageDiff数据框。head()函数返回exponentialMovingAverageDiff数据框的前部分。exponentialMovingAverageDiff数据框作为输入参数传递：

> head(exponentialMovingAverageDiff)

结果如下：

探索价格变化的总结：为此使用summary()函数。此函数提供一系列描述性统计，以生成exponentialMovingAverageDiff数据框的结果摘要。

> summary(exponentialMovingAverageDiff)

结果如下：

我们现在将比较GSPC系列快速移动平均与GSPC系列的慢速移动平均。为此，将GSPC作为价格矩阵传递。fast = 12代表快速移动平均的周期数，slow = 26代表慢速移动平均的周期数，signal = 9代表移动平均的信号：

> MACD <- MACD(Op(GSPC),fast = 12, slow = 26, signal = 9)

打印 MACD 数据框：tail() 函数返回 MACD 数据框的最后部分。MACD 数据框作为输入参数传递：

> tail(MACD)

结果如下：

使用 summary() 函数探索价格变化摘要：

> summary(MACD)

结果如下：

接下来，我们将抓取信号线作为指标。结果存储在 MACDsignal 数据框中：

> MACDsignal <- MACD[,2]

计算 布林带：它们是范围指标，从移动平均线计算标准差。布林带在以下逻辑下运行：货币对的价格最有可能趋向于其平均值；因此，当它偏离太多，比如说两个标准差之外时，它将回落到其移动平均线。BBands() 函数用于计算布林带。GSPC 作为对象传递，n=20 表示移动平均期的数量。sd=2 表示两个标准差：

> BollingerBands <- BBands(Op(GSPC),n=20,sd=2)

现在让我们打印 BollingerBands 数据框：

> tail(BollingerBands)

结果如下：

探索价格变化的摘要：

> summary(BollingerBands)

结果如下：

现在让我们从 BollingerBands 抓取信号线作为指标：

> PercentageChngpctB <- BollingerBands[,4]

打印 PercentageChngpctB 数据框：

> tail(PercentageChngpctB)

结果如下：

探索 PercentageChngpctB 的变化摘要：

> summary(PercentageChngpctB)

结果如下：

查找收盘价和开盘价之间的差异：

> Price <- Cl(GSPC)-Op(GSPC)

打印 price 数据框：

> tail(Price)

结果如下：

结合 relativeStrengthIndex3、expMvAvg5Cross、MACDsignal、PercentageChngpctB 和 Price 数据框：结果随后存储在 DataSet 数据框中：

> DataSet <- data.frame(relativeStrengthIndex3, expMvAvg5Cross, MACDsignal, PercentageChngpctB, Price)

探索 DataSet 数据框的内部结构：str() 函数显示数据框的内部结构。DataSet 作为 R 对象传递给 str() 函数：

> str(DataSet)

结果如下：

计算指标、创建数据集和移除点：

> DataSet <- DataSet[-c(1:33),]

探索 DataSet 数据框的维度：dim() 函数返回 DataSet 数据框的维度。DataSet 数据框作为输入参数传递。结果清楚地表明有 1,176 行数据和 5 列：

> dim(DataSet)

结果如下：

命名列：c() 函数用于将参数组合成向量：

> colnames(DataSet) <- c("RSI3","EMAcross","MACDsignal","BollingerB","Price")

探索 DataSet 数据框的维度：

> str(DataSet)

结果如下：

步骤 4 - 为模型构建准备数据

将数据集归一化到 0 到 1 之间：

> Normalized <- function(x) {(x-min(x))/(max(x)-min(x))}

调用归一化数据集的函数：

> NormalizedData <- as.data.frame(lapply(DataSet,Normalized))

打印 NormalizedData 数据框：

> tail(NormalizedData)

结果如下：

构建训练数据集：NormalizedData 数据框中的 1:816 数据元素将被用作训练数据集。训练数据集应存储在 TrainingSet 中：

> TrainingSet <- NormalizedData[1:816,]

探索 TrainingSet 数据框的维度：

> dim(TrainingSet)

结果如下：

探索 TrainingSet 的变化摘要：

> summary(TrainingSet)

结果如下：

构建测试数据集：NormalizedData 数据框中的 817:1225 数据元素将被用作训练数据集。此测试数据集应存储在 TestSet 中：

> TestSet <- NormalizedData[817:1225 ,]

探索 TrainingSet 数据框的维度：

> dim(TestSet)

结果如下：

探索 TestSet 的变化摘要：

> summary(TestSet)

结果如下：

第 5 步 - 构建模型

构建神经网络：neuralnet() 函数使用不带权重回溯的反向传播算法训练神经网络。Price~RSI3+EMAcross+MACDsignal+BollingerB 是要拟合的模型的描述。data=TrainingSet 是包含公式中指定变量的数据框。hidden=c(3,3) 指定了每层的隐藏神经元（顶点）数量。learningrate=.001 表示反向传播算法使用的学习率。algorithm="backprop" 指的是反向传播算法：

> nn1 <- neuralnet(Price~RSI3+EMAcross+MACDsignal+BollingerB,data=TrainingSet, hidden=c(3,3), learningrate=.001,algorithm="backprop")

绘制神经网络：

> plot(nn1)

结果如下：

测量失业率

失业率定义为失业的劳动力占总劳动力的百分比，且积极寻找工作并愿意工作。根据国际劳工组织（ILO）的定义，失业者是指积极寻找工作但没有工作的人。失业率是衡量失业人数同时失业和寻找工作的指标。

准备就绪

为了使用神经网络进行失业率测量，我们将使用收集到的威斯康星州失业率数据集。

第 1 步 - 收集和描述数据

为此，我们将使用标题为 FRED-WIUR.csv 的 CSV 数据集。有 448 行数据。有两个数值变量如下：

• 日期

• 值

此数据集显示了 1976 年 1 月 1 日至 2013 年 4 月 1 日间威斯康星州的失业率。

如何操作...

让我们深入了解。

第 2 步 - 探索数据

首先，需要加载以下包：

 > install.packages("forecast ")
 > install.packages("lmtest") 
 > install.packages("caret ")
 > library(forecast)
 > library(lmtest)
 > library(caret)

注意

版本信息：本页面的代码在 R 版本 3.3.0 中进行了测试

让我们探索数据并了解变量之间的关系。我们将首先导入名为 FRED-WIUR.csv 的 CSV 数据文件。我们将数据保存到 ud 数据框中：

> ud <- read.csv("d:/FRED-WIUR.csv", colClasses=c('Date'='Date'))

打印 ud 数据框：tail() 函数返回 ud 数据框的最后部分。将 ud 数据框作为输入参数传递：

> tail(ud)

结果如下：

命名列：使用 c() 函数将参数组合成向量：

> colnames(ud) <- c('date', 'rate')

使用 as.Date() 函数将字符表示和 Date 类的对象转换为日期：

> ud$date <- as.Date(ud$date)

探索失业数据的摘要：为此，使用 summary() 函数。该函数提供一系列描述性统计，以生成 ud 数据框的结果摘要：

> summary (ud)

结果如下：

现在，让我们从第 1 行到第 436 行创建基础数据：

> ud.b <- ud[1:436,]

探索基础失业数据的摘要。为此，使用 summary() 函数。该函数提供一系列描述性统计，以生成 ud.b 数据框的结果摘要：

> summary(ud.b)

结果如下：

现在，让我们从第 437 行到第 448 行创建测试数据：

> ud.p <- ud[437:448,]

探索测试失业数据的摘要：

> summary(ud.p)

结果如下：

从 1976 年创建基础时间序列数据：ts() 函数作为创建时间序列对象的函数。ud.b$rate 代表观察到的时序值向量：

> ud.ts <- ts(ud.b$rate, start=c(1976, 1), frequency=12)

打印 ud.ts 数据框的值：

> ud.ts

结果如下：

创建测试时间序列数据：ts() 函数创建时间序列对象。ud.b$rate 代表观察到的时序值向量：

> ud.p.ts <- ts(ud.p$rate, start=c(2012, 5), frequency=12)

打印 ud.ts 数据框的值：

> ud.ts

结果如下：

绘制基础时间序列数据：

> plot.ts(ud.ts)

结果如下：

绘制测试时间序列数据：

> plot.ts(ud.p.ts)

结果如下：

第三步 - 准备和验证模型

计算基础时间序列数据集的平均值。meanf() 函数返回对 ud.ts 数据集应用 i.i.d 模型后的预测和预测区间。12 表示预测的周期：

> mean <- meanf(ud.ts, 12)

对具有漂移的基础时间序列进行预测和预测区间。rwf() 函数对时间序列 ud.ts 进行随机游走预测并返回。参数 12 表示预测的周期：

> forecast_randomwalk <- rwf(ud.ts, 12)

从 ARIMA(0,0,0)(0,1,0)m 基础时间序列对随机游走进行预测和预测区间：snaive()函数对时间序列ud.ts执行 ARIMA(0,0,0)(0,1,0)m，并返回预测结果。参数12表示预测的周期：

> forecast_arima <- snaive(ud.ts, 12)

预测基础时间序列的漂移。rwf()函数对时间序列ud.ts上的随机游走进行预测并返回结果。参数12表示预测的周期。drift=T是一个逻辑标志，用于拟合带有漂移模型的随机游走：

> drift <- rwf(ud.ts, 12, drift=T)

接下来，我们将为趋势基础时间序列数据准备线性拟合模型。tslm()函数将线性模型拟合到ud.ts时间序列。ud.ts~trend公式表示必须考虑趋势成分：

> m1 <- tslm(ud.ts~trend)

为基础时间序列数据准备趋势和季节性的线性拟合模型：tslm()函数将线性模型拟合到ud.ts时间序列。ud.ts~trend+season公式表示趋势和季节性成分必须被考虑：

> m2 <- tslm(ud.ts~trend+season)

residuals()是一个通用函数，在为趋势基础时间序列数据拟合模型后，从对象m1中提取模型残差。

> residual_1 <- residuals(m1)

绘制残差模型：

> plot(residual_1, ylab="Residuals",xlab="Year", title("Residual - Trends"), col = "red")

结果如下：

现在我们来看如何估计自协方差函数。residual_1是单变量数值时间序列对象：

> acf(residual_1, main="ACF of residuals")

结果如下：

residuals()是一个通用函数，在为趋势基础时间序列数据拟合模型后，从对象m2中提取模型残差。

> residual_2 <- residuals(m2)

绘制残差模型：

> plot(residual_2, ylab="Residuals",xlab="Year",title("Residual - Trends + Seasonality"), col = "red")

结果如下：

> acf(residual_2, main="ACF of residuals")

结果如下：

杜宾-沃森检验用于确定线性回归或多元回归的残差是否独立。在杜宾-沃森检验中通常考虑的假设如下：

测试统计量如下：

在此方程中，，是单个的观测值，而是单个的预测值。

随着序列相关性的增加，的值降低。对于的不同值（解释变量的数量）和，已经为上、下临界值和编制了表格：

如果拒绝

如果不拒绝

如果的测试结果不确定。

对基础时间序列数据的趋势进行线性拟合模型的 Durbin-Watson 测试：

> dwtest(m1, alt="two.sided")

结果如下：

对基础时间序列数据的趋势和季节性进行线性拟合模型的 Durbin-Watson 测试：

 > dwtest(m2, alt="two.sided")

结果如下：

使用 LOESS 将基础数据时间序列分解为周期、季节、趋势和不规则成分：

> m3 <- stl(ud.ts, s.window='periodic')

绘制分解后的基础数据时间序列图：

> plot(m3)

结果如下：

对基础数据时间序列执行指数平滑状态空间模型。ets()函数返回ud.ts时间序列上的ets模型。ZZZ - "Z"表示自动选择。第一个字母表示误差类型，第二个字母表示趋势类型，第三个字母表示季节类型：

> m4 <- ets(ud.ts, model='ZZZ')

绘制基础数据时间序列的指数平滑状态空间模型图：

> plot(m4)

结果如下：

返回基础数据时间序列单变量 ARIMA 的阶数：

> m5 <- auto.arima(ud.ts)

绘制基础数据时间序列的单变量 ARIMA 图：

> plot(forecast(m5, h=12))

结果如下：

构建前馈神经网络模型：nnetar()函数使用单个隐藏层和滞后输入构建前馈神经网络，用于预测基础数据的单变量时间序列：

> m6 <- nnetar(ud.ts)

打印前馈神经网络模型的值：

> m6

结果如下：

绘制前馈神经网络模型图：

> plot(forecast(m6, h=12))

结果如下：

步骤 4 - 预测和测试构建的模型的准确性

使用测试数据时间序列测试基础数据时间序列平均值的准确性。accuracy()函数返回预测准确性的汇总度量范围。ud.p.ts是测试数据时间序列：

> a1 <- accuracy(mean, ud.p.ts)

使用漂移测试预测和预测的基础数据时间序列的准确性：

> a2 <- accuracy(forecast_randomwalk, ud.p.ts)

使用 ARIMA(0,0,0)(0,1,0)m 测试预测和预测的基础数据时间序列的准确性：

> a3 <- accuracy(forecast_arima, ud.p.ts)

测试基础数据时间序列漂移的准确性：

> a4 <- accuracy(drift, ud.p.ts)

将结果组合到表格中：

> a.table <- rbind(a1, a2, a3, a4)

打印结果：

> a.table

结果如下：

预测趋势的基础时间序列数据的线性拟合模型。h=12表示预测的周期：

> f1 <- forecast(m1, h=12)

预测趋势和季节性的基础时间序列数据的线性拟合模型：

> f2 <- forecast(m2, h=12)

使用 LOESS 将分解的基础数据时间序列预测为周期、季节、趋势和不规则成分：

> f3 <- forecast(m3, h=12)

预测基础数据时间序列的指数平滑状态空间模型：

> f4 <- forecast(m4, h=12)

预测基础数据时间序列的有序单变量 ARIMA：

> f5 <- forecast(m5, h=12)

预测具有单个隐藏层的前馈神经网络模型：

> f6 <- forecast(m6, h=12)

测试预测的趋势的基础时间序列数据的线性拟合模型的准确性：

> a5 <- accuracy(f1, ud.p.ts)

测试预测的趋势和季节性的基础时间序列数据的线性拟合模型的准确性：

> a6 <- accuracy(f2, ud.p.ts)

使用 LOESS 测试预测分解的基础数据时间序列（周期、季节、趋势和不规则成分）的准确性：

> a7 <- accuracy(f3, ud.p.ts)

测试基础数据时间序列预测的指数平滑状态空间模型的准确性：

> a8 <- accuracy(f4, ud.p.ts)

测试基础数据时间序列预测的有序单变量 ARIMA 的准确性：

> a9 <- accuracy(f5, ud.p.ts)

测试预测的前馈神经网络模型（具有单个隐藏层）的准确性：

> a10 <- accuracy(f6, ud.p.ts)

将结果组合到表格中：

> a.table.1 <- rbind(a5, a6, a7, a8, a9, a10)

打印结果：

> a.table.1

结果如下：

第十一章：深度学习

在本章中，我们将介绍以下内容：

循环神经网络 - 预测周期信号

简介

大多数机器学习算法由于预定义的表示和输入特征而表现良好。机器学习算法通过优化权重来最佳地做出最终预测，而表示学习试图自动学习良好的特征或表示。深度学习算法通过增加复杂性来尝试在多个表示级别上学习。深度架构由多个非线性行为级别组成，例如具有许多隐藏层的神经网络。深度学习技术的目标主要是学习特征层次。深度学习技术可以分为三大类；用于无监督或生成学习的深度网络，用于监督学习的深度网络和混合深度网络

循环神经网络 - 预测周期信号

振荡器是产生特定、周期性波形的电路，如方波、三角波、锯齿波和正弦波。为了生成输出，振荡器通常使用某种形式的主动设备-灯，它被电阻器、电容器和电感器所包围。振荡器的两种主要类型是弛豫振荡器和正弦振荡器。三角波、锯齿波和其他非正弦波形是通过弛豫振荡器产生的，而正弦振荡器由外部组件和放大器组成，以产生振荡。通常，纯正弦波中不含有谐波，它们只包含一个频率。

准备工作...

任务是从一个带噪声的正弦波中预测余弦波。使用 5Hz 频率的正弦波，其中包含一些正态分布的噪声和一个平滑的余弦波。创建的数据集是一组 10 个序列，每个序列包含 40 个观测值。

如何做...

需要在第一步加载以下包：

 > install.packages("rnn")
 > library(rnn)

将初始种子设置为随机数，以实现可重复性：

> set.seed(10)

初始化所需的频率：

> f <- 5

创建所需的向量：

> w <- 2*pi*f

生成序列：seq()函数生成常规序列。0.005是起始值，2是结束值。by=0.005确定增量序列：

> t <- seq(0.005,2,by=0.005)

生成sin和cos值：

 > x <- sin(t*w) + rnorm(200, 0, 0.25)
 > y <- cos(t*w)

生成时间序列样本：matrix()函数从x和y值创建矩阵。nrow = 40表示所需的行数：

 > X <- matrix(x, nrow = 40)
 > Y <- matrix(y, nrow = 40)

绘制带噪声的波形：plot()函数是用于绘制 R 对象的通用函数。as.vector(X)数据框作为函数值传递。type='l'表示线条：

> plot(as.vector(X), col='blue', type='l', ylab = "x-matrix, y-matrix", main = "Noisy waves")

结果如下：

> lines(as.vector(Y), col = "red")

结果如下：

标准化X的值。值的范围在 0 到 1 之间：

> X <- (X - min(X)) / (max(X) - min(X))

打印X的值：

> X

结果如下：

标准化 Y 的值。值的范围介于 0 和 1 之间：

> X <- (X - min(X)) / (max(X) - min(X))

打印 X 的值：

> X

结果如下：

转置 X 和 Y 的值：

 > X <- t(X)
 > Y <- t(Y)

创建训练集和测试集：

 > train <- 1:8
 > test <- 9:10

训练循环神经网络。Y = Y[train,] 表示输出值的数组。X = X[train,] 表示输入值的数组。learningrate = 0.05 表示权重迭代的速率。hidden_dim = 16 是隐藏层的维度。numepochs = 1500 是整个数据集进行训练的次数。

这个阶段将花费时间。所需时间取决于学习率、维度数量以及整个数据集进行训练的次数：

> model <- trainr(Y = Y[train,],X = X[train,],learningrate = 0.05,hidden_dim = 16,numepochs = 1500)

结果如下：

预测循环神经网络（Recurrent Neural Network）的输出：

> Y_predicted <- predictr(model, X)

绘制 实际值与预测值 的对比图。输出构成训练集和测试集：

> plot(as.vector(t(Y)), col = 'red', type = 'l', main = "Actual values vs Predicted values", ylab = "Y, Y-predicted")

结果如下：

> lines(as.vector(t(Y_predicted)), type = 'l', col = 'blue')

结果如下：

绘制 实际值与预测值 的对比图。输出仅构成测试集：

> plot(as.vector(t(Y[test,])), col = 'red', type='l', main = "Actual vs predicted: testing set", ylab = "Y,Y-predicted")

结果如下：

> lines(as.vector(t(Y_predicted[test,])), type = 'l', col = 'blue')

结果如下：

第十二章. 案例研究 - 探索世界银行数据

简介

世界银行指标（WDI）是世界银行汇编的关于全球发展增长及其对人们生活质量的国际可比和可衡量统计数据的集合。通过分析来自 200 多个经济体和 50 多个合作伙伴组织收集的数据，衡量各国、地区和收入群体的发展状况，展示了 1400 多个指标。2015 年 9 月 25 日，联合国大会正式通过了 2030 年可持续发展议程，以指导未来 15 年的全球行动。可持续发展目标（SDG）的五大主要关注主题是人、地球、繁荣、和平和伙伴关系。各国已承诺消除贫困和饥饿，并确保所有人都能在尊严和平等的环境中充分发挥其潜能；保护地球免受退化，并就气候变化采取紧急行动；确保所有人都能享有繁荣充实的生活，并且进步与自然和谐共生；培育和平、公正和包容的社会，消除恐惧和暴力；并动员实施 2030 年议程的手段，重点关注最贫困和最脆弱的人群，通过强大的全球伙伴关系。对于这 17 个目标，世界银行发展数据组、全球实践和跨领域解决方案领域的专家为每个目标选择了指标，以识别和分析重要趋势和挑战，并就测量问题进行讨论。《世界发展指标》是众多国际机构、200 多个国家统计办公室以及许多其他机构合作的结果。

探索世界银行数据

2012 年，全球 13%的人口生活在每天 1.90 美元的国际贫困线以下，而 1990 年这一比例为 37%。所有地区的下降都为提前实现全球减半极端贫困的千年发展目标做出了贡献。目标是到 2030 年消除所有形式的贫困，并为贫困人口提供社会保障，增加基本服务的获取，并支持受冲突和气候相关灾害影响的人。

在低收入国家记录的死亡中，超过一半是由于传染病或母体、围产期或营养状况。而在中高收入国家，记录的死亡中超过三分之二是由于非传染性疾病。全球能源使用量在 1990 年至 2013 年之间增长了约 54%。能源获取是发展的基础，但随着经济的演变，收入增长和人口增长对能源的需求也在增加。能源，尤其是电力，对于提高低收入和中收入国家人们的生活水平至关重要。

准备中...

为了执行再保险合同的定价，我们将使用从飓风数据集收集的数据集。

为了对世界银行的数据模式进行分析，我们将使用从以下数据集收集的数据集：

• 全球总人口（1960-2015）

• 所有国家和地区的出生预期寿命（1960-2014）

• 所有国家和地区的生育率（每名女性的出生数），时间范围为（1960-2014）

• 所有国家和地区的 GDP（以美元计），时间范围为（1960-2015）

• 所有国家和地区的贫困人口比例（1960-2016）

• 所有国家和地区的卫生设施普及率（1960-2016）

• 所有国家和地区的电力普及率（1960-2016）

• 所有国家和地区的二氧化碳排放量（1960-2016）

第 1 步 - 收集和描述数据

世界银行用于分析的数据集可以从世界银行数据库免费下载。

如何操作...

让我们深入了解细节。

第 2 步 - 下载数据

加载以下包：

 > install.packages("wbstats")
 > install.packages("data.table")
 > install.packages("googleVis")

备注

版本信息：本页面的代码在 R 版本 3.3.0（2016-05-03）上进行了测试

以下每个库都需要安装：

 > library(wbstats)
 > library(data.table)
 > library(googleVis)

让我们下载数据并了解变量之间的关系。我们将首先从世界银行网站下载数据。data.table()函数允许快速聚合大型数据集，排序连接，通过组添加/修改/删除列，列出列，友好的文件读取器，以及并行文件写入器。wb()函数使用世界银行 API 下载所需信息。indicator 代表指标代码的字符向量。

指标代码如下：

• SP.POP.TOTL: 全球总人口（1960-2015）

• SP.DYN.LE00.IN: 所有国家和地区的出生预期寿命（1960-2014）

• SP.DYN.TFRT.IN: 所有国家和地区的生育率（每名女性的出生数），时间范围为（1960-2014）

结果存储在Pop_LifeExp_FertRt数据框中。使用以下命令：

> Pop_LifeExp_FertRt <- data.table(wb(indicator = c("SP.POP.TOTL", "SP.DYN.LE00.IN", "SP.DYN.TFRT.IN"), startdate = 1960, enddate = 2016))

指标代码如下：

• SP.POP.TOTL: 全球总人口（1960-2015）

• NY.GDP.MKTP.CD-GDP: 所有国家和地区的 GDP（以美元计），时间范围为（1960-2015）

• SI.POV.2DAY: 所有国家和地区的贫困人口比例（1960-2016）

结果存储在Pop_GDPUSD_HeadCnt数据框中。使用以下命令：

> Pop_GDPUSD_HeadCnt <- data.table(wb(indicator = c("SP.POP.TOTL", "NY.GDP.MKTP.CD", "SI.POV.2DAY"), startdate = 1960, enddate = 2016))

指标代码如下：

• SP.POP.TOTL: 全球总人口（1960-2015）

• NY.GDP.MKTP.CD: 所有国家和地区的 GDP（以美元计），时间范围为（1960-2015）

• SH.STA.ACSN: 所有国家和地区的卫生设施普及率（1960-2016）

结果存储在Pop_GDPUSD_Sanitation数据框中。使用以下命令：

> Pop_GDPUSD_Sanitation <- data.table(wb(indicator = c("SP.POP.TOTL", "NY.GDP.MKTP.CD", "SH.STA.ACSN"), startdate = 1960, enddate = 2016))

指标代码如下：

• NY.GDP.MKTP.CD: 所有国家和地区的 GDP（以美元计），时间范围为（1960-2015）

• EG.ELC.ACCS.ZS: 所有国家和地区的电力普及率（1960-2016）

• EN.ATM.CO2E.KT：所有国家和地区的每人大电力消耗 KWh（1960-2016）

结果存储在 GDPUSD_Electricity_CO2 数据框中。使用以下命令：

> GDPUSD_Electricity_CO2 <- data.table(wb(indicator = c("NY.GDP.MKTP.CD", "EG.ELC.ACCS.ZS", "EN.ATM.CO2E.KT"), startdate = 1960, enddate = 2016))

步骤 3 - 探索数据

探索 Pop_LifeExp_FertRt 数据框的维度：dim() 函数返回 Pop_LifeExp_FertRt 数据框的维度。Pop_LifeExp_FertRt 数据框作为输入参数传递。结果明确指出有 41150 行数据和六个列：

> dim(Pop_LifeExp_FertRt)

结果如下：

探索 Pop_LifeExp_FertRt 数据框的维度：结果明确指出有 27023 行数据和六个列：

> dim(Pop_GDPUSD_HeadCnt)

结果如下：

探索 Pop_GDPUSD_Sanitation 数据框的维度：结果明确指出有 31884 行数据和六个列：

> dim(Pop_GDPUSD_Sanitation)

结果如下：

探索 GDPUSD_Electricity_CO2 数据框的维度：结果明确指出有 23994 行数据和六个列：

> dim(GDPUSD_Electricity_CO2)

结果如下：

探索 Pop_LifeExp_FertRt 数据框的内部结构：str() 函数显示数据框的内部结构。Pop_LifeExp_FertRt 作为 R 对象传递给 str() 函数：

> str(Pop_LifeExp_FertRt)

结果如下：

探索 Pop_GDPUSD_HeadCnt 数据框的内部结构：

> str(Pop_GDPUSD_HeadCnt)

结果如下：

探索 Pop_GDPUSD_Sanitation 数据框的内部结构：

> str(Pop_GDPUSD_Sanitation)

结果如下：

探索 GDPUSD_Electricity_CO2 数据框的内部结构：

> str(GDPUSD_Electricity_CO2)

结果如下：

探索 GDPUSD_Electricity_CO2 数据框的内部结构：

> str(GDPUSD_Electricity_CO2)

结果如下：

打印 Pop_LifeExp_FertRt 数据框：head() 函数返回 Pop_LifeExp_FertRt 数据框的前部分。Pop_LifeExp_FertRt 数据框作为输入参数传递：

> head(Pop_LifeExp_FertRt)

结果如下：

打印 Pop_GDPUSD_HeadCnt 数据框：

> head(Pop_GDPUSD_HeadCnt)

结果如下：

打印 Pop_GDPUSD_Sanitation 数据框：

> head(Pop_GDPUSD_Sanitation)

结果如下：

打印 GDPUSD_Electricity_CO2 数据框：

> head(GDPUSD_Electricity_CO2)

结果如下：

探索 SP.POP.TOTL 数据框的维度：dim() 函数返回 SP.POP.TOTL 数据框的维度。SP.POP.TOTL 数据框作为输入参数传递。结果清楚地表明有 14623 行数据和六个列：

> dim(wb(indicator = "SP.POP.TOTL"))

结果如下：

探索 SP.DYN.LE00.IN 数据框的维度：

> dim(wb(indicator = "SP.DYN.LE00.IN"))

结果如下：

探索 SP.DYN.TFRT.IN 数据框的维度：

> dim(wb(indicator = " SP.DYN.TFRT.IN "))

结果如下：

探索 NY.GDP.MKTP.CD 数据框的维度：

> dim(wb(indicator = " NY.GDP.MKTP.CD"))

结果如下：

探索 SI.POV.2DAY 数据框的维度：

> dim(wb(indicator = " SI.POV.2DAY "))

结果如下：

探索 SH.STA.ACSN 数据框的维度：

> dim(wb(indicator = " SH.STA.ACSN "))

结果如下：

探索 EG.ELC.ACCS.ZS 数据框的维度：

> dim(wb(indicator = "EG.ELC.ACCS.ZS"))

结果如下：

探索 EN.ATM.CO2E.KT 数据框的维度：

> dim(wb(indicator = "EN.ATM.CO2E.KT"))

结果如下：

使用 wbcountries() 函数从世界银行 API 下载更新的国家和区域信息：

> Countries <- data.table(wbcountries())

打印 Countries 数据框：head() 函数返回 Countries 数据框的前部分：

> head(Countries)

结果如下：

步骤 4 - 构建模型

对 Pop_LifeExp_FertRt 数据表进行排序：setkey() 函数对 Pop_LifeExp_FertRt 数据表进行排序并标记为已排序。排序的列是键。键位于 iso2c 列；iso2c 列始终按升序排序。表通过引用进行更改，因此非常节省内存：

> setkey(Pop_LifeExp_FertRt, iso2c)

对 Pop_GDPUSD_HeadCnt 数据表进行排序：

> setkey(Pop_GDPUSD_HeadCnt, iso2c)

对 Pop_GDPUSD_Sanitation 数据表进行排序：

> setkey(Pop_GDPUSD_Sanitation, iso2c)

对 GDPUSD_Electricity_CO2 数据表进行排序：

> setkey(GDPUSD_Electricity_CO2, iso2c)

对 Countries 数据表进行排序：

> setkey(Countries, iso2c)

打印 Countries 数据表：

> head(setkey(Countries, iso2c))

结果如下：

在数据集中添加区域的同时从 Pop_LifeExp_FertRt 数据集中移除聚合：

> Pop_LifeExp_FertRt <- Countries[Pop_LifeExp_FertRt][ ! region %in% "Aggregates"]

打印 Pop_LifeExp_FertRt 数据表：

> head(Pop_LifeExp_FertRt)

结果如下：

在数据集中添加区域的同时从 Pop_GDPUSD_HeadCnt 数据集中移除聚合：

> Pop_GDPUSD_HeadCnt <- Countries[Pop_GDPUSD_HeadCnt][ ! region %in% "Aggregates"]

在数据集中添加区域的同时从 Pop_GDPUSD_Sanitation 数据集中移除聚合：

> Pop_GDPUSD_Sanitation <- Countries[Pop_GDPUSD_Sanitation][ ! region %in% "Aggregates"]

在数据集中添加区域的同时从 GDPUSD_Electricity_CO2 数据集中移除聚合：

> GDPUSD_Electricity_CO2 <- Countries[GDPUSD_Electricity_CO2][ ! region %in% "Aggregates"]

 > wPop_LifeExp_FertRt <- reshape(Pop_LifeExp_FertRt[, list(country, region, date, value, indicator)], v.names = "value", idvar=c("date", "country", "region"), timevar="indicator", direction = "wide")
 > wPop_GDPUSD_HeadCnt <- reshape(Pop_GDPUSD_HeadCnt[, list(country, region, date, value, indicator)], v.names = "value", idvar=c("date", "country", "region"), timevar="indicator", direction = "wide")
 > wPop_GDPUSD_Sanitation <- reshape(Pop_GDPUSD_Sanitation[, list(country, region, date, value, indicator)], v.names = "value", idvar=c("date", "country", "region"), timevar="indicator", direction = "wide")
 > wGDPUSD_Electricity_CO2 <- reshape(GDPUSD_Electricity_CO2[, list(country, region, date, value, indicator)], v.names = "value", idvar=c("date", "country", "region"), timevar="indicator", direction = "wide")

打印数据框 wPop_LifeExp_FertRt 的内容：

> wPop_LifeExp_FertRt

结果如下：

打印数据框 wGDPUSD_Electricity_CO2 的内容：

> wGDPUSD_Electricity_CO2

结果如下：

将wPop_LifeExp_FertRt、wPop_GDPUSD_HeadCnt、wPop_GDPUSD_Sanitation和wGDPUSD_Electricity_CO2数据集从字符格式转换为整数格式：

 > wPop_LifeExp_FertRt[, date := as.integer(date)]
 > wPop_GDPUSD_HeadCnt[, date := as.integer(date)]
 > wPop_GDPUSD_Sanitation[, date := as.integer(date)]
 > wGDPUSD_Electricity_CO2[, date := as.integer(date)]

设置名称：setnames()函数设置wPop_LifeExp_FertRt、wPop_GDPUSD_HeadCnt、wPop_GDPUSD_Sanitation和wGDPUSD_Electricity_CO2对象的名字：

    > setnames(wPop_LifeExp_FertRt, names(wPop_LifeExp_FertRt), c("Country", "Region", "Year", "Population", "Fertility", "LifeExpectancy"))
    > setnames(wPop_GDPUSD_HeadCnt, names(wPop_GDPUSD_HeadCnt), c("Country", "Region", "Year", "Population", "GDPUSD", "PovertyHead"))
    > setnames(wPop_GDPUSD_Sanitation, names(wPop_GDPUSD_Sanitation), c("Country", "Region", "Year", "Population", "GDPUSD", "SanitationAccess"))
    > setnames(wGDPUSD_Electricity_CO2, names(wGDPUSD_Electricity_CO2), c("Country", "Region", "Year", "GDPUSD", "ElectricityConsumption", "CO2Emissions"))

第 5 步 - 绘制模型

按以下步骤绘制wPop_LifeExp_FertRt数据框模型。gvisMotionChart()函数读取wPop_LifeExp_FertRt数据框。它使用 Google Visualization API 创建包含在网页中的文本输出。图表由网络浏览器在 Flash 中渲染。动态的动图探索指标。wPop_LifeExp_FertRt是数据框。idvar = "Country"表示要分析的数据的列名。timevar = "Year"是表示时间维度的数据的列名。xvar = "LifeExpectancy"是要绘制在x轴上的数据的数值向量。yvar = "Fertility"是要绘制在 y 轴上的数据的数值向量。sizevar = "Population"表示要映射到实际像素值的列值。colorvar = "Region"标识气泡。使用以下命令：

 > pltPop_LifeExp_FertRt <- gvisMotionChart(wPop_LifeExp_FertRt, idvar = "Country", timevar = "Year", xvar = "LifeExpectancy", yvar = "Fertility", sizevar = "Population", colorvar = "Region")
 > plot(pltPop_LifeExp_FertRt)

绘制wPop_GDPUSD_HeadCnt数据框模型：

 > pltPop_GDPUSD_HeadCnt <- gvisMotionChart(wPop_GDPUSD_HeadCnt, idvar = "Country", timevar = "Year", xvar = "GDPUSD", yvar = "PovertyHead", sizevar = "Population", colorvar = "Region")
 > plot(pltPop_GDPUSD_HeadCnt)

绘制wPop_GDPUSD_Sanitation数据框模型：

 > pltPop_GDPUSD_Sanitation <- gvisMotionChart(wPop_GDPUSD_Sanitation, idvar = "Country", timevar = "Year", xvar = "GDPUSD", yvar = "SanitationAccess", sizevar = "Population", colorvar = "Region")
 > plot(pltPop_GDPUSD_Sanitation)

绘制pltGDPUSD_Electricity_CO2数据框模型：

 > pltGDPUSD_Electricity_CO2 <- gvisMotionChart(wGDPUSD_Electricity_CO2, idvar = "Country", timevar = "Year", xvar = "GDPUSD", yvar = "ElectricityAccess", sizevar = "CO2Emissions", colorvar = "Region")
 > plot(pltGDPUSD_Electricity_CO2)

结果如下，生育率与预期寿命的关系：

人口增长：

以美元计算的 GDP 增长：

贫困人口比例与人口增长的关系：

人群获得卫生设施的增长：

卫生设施普及率：

卫生设施普及率的改善与人口增长的关系：

所有国家和地区的电力普及率人口：

二氧化碳排放（对数刻度）：

二氧化碳排放与电力消费的关系：

第十三章. 案例研究 - 重新保险合同的定价

简介

如同其名，再保险是从保险业务发展而来的，其使用的程度不仅取决于直接承保人要承保的风险的金额，还取决于这些风险的特征。可以进行的再保险业务的量主要取决于任何给定时间可用的直接业务量。再保险的想法根植于与保险产生相同的人类本能，即希望一个人的损失由许多人分担。

重新保险合同的定价

保险公司安排重新保险的关键目标包括但不限于：通过将不会因财务限制而承担的那部分风险转嫁给再保险公司，从而增加处理更大风险的能力；增强接受比资本允许的更大保额的能力；通过再保险公司吸收更大的索赔或灾难损失，使年度经营结果稳定；通过加强承保人建立既在规模又在风险质量上同质化的账户的尝试，增加盈利的机会；能够承保新的风险敞口。再保险的功能可以视为提供保护增加能力、财务稳定性、稳定索赔比率、积累不同类别的索赔、风险分散、保护偿付能力边际和稳定利润的服务。再保险有助于吸收由经济变化、社会变化、保险方法的变化以及科学发展引起的新的风险敞口。重新保险合同可以安排的方式只有两种：一种是自愿再保险，针对单一保单的一次性；另一种是合同再保险，针对定义好的保单组的自动安排。

准备中...

为了进行重新保险合同的定价，我们将使用收集在飓风数据集上的数据集。

第 1 步 - 收集和描述数据

将使用名为 publicdatamay2007.xls 的 XLS 格式数据集。该数据集采用标准格式。共有 207 行数据。有 7 个变量。数值变量如下：

• 年份

• 基础经济损失

• 标准化 PL05

• 标准化 CL05

非数值变量如下：

• 飓风描述

• 状态

• 类别

如何操作...

让我们深入了解细节。

第 2 步 - 探索数据

加载以下包：

 > install.packages("gdata")
 > install.packages("evir")
 > library(gdata)
 > library(evir)

注意

版本信息：本页面的代码在 R 版本 3.2.2 中进行了测试

让我们探索数据并了解变量之间的关系，如下所示。我们将首先导入名为 publicdatamay2007.xls 的 XLS 数据文件。我们将把数据保存到 StormDamageData 数据框中：

> StormDamageData <- read.xls("d:/publicdatamay2007.xls", sheet = 1)

打印 StormDamageData 数据框：head() 函数返回 StormDamageData 数据框的第一部分。将 StormDamageData 数据框作为输入参数传递：

> head(StormDamageData)

结果如下：

tail()函数返回StormDamageData框的最后部分，如下所示。将StormDamageData框作为输入参数传递。

> tail(StormDamageData)

结果如下：

探索StormDamageData数据框的维度：dim()函数返回StormDamageData框的维度。将StormDamageData数据框作为输入参数传递。结果清楚地表明有 207 行数据和九列：

> dim(StormDamageData)

结果如下：

步骤 3 - 计算个别损失索赔

格式化数据：包装函数ChangeFormat消除了传递的值中的逗号（,），并以数值形式返回结果：

 > ChangeFormat <- function(x){
 x = as.character(x)
 for(i in 1:10){x=sub(",","",as.character(x))}
 return(as.numeric(x)) }

将StormDamageData数据框存储在基础中：

> base <- StormDamageData[,1:4]

调用包装函数，ChangeFormat：将StormDamageData数据框中的Base.Economic.Damage作为输入传递。Vectorize()函数创建了对ChangeFormat()函数的包装。结果存储在base$Base.Economic.Damage数据框中：

> base$Base.Economic.Damage <- Vectorize(ChangeFormat)(StormDamageData$Base.Economic.Damage)

调用包装函数，ChangeFormat：将StormDamageData数据框中的Normalized.PL05作为输入传递。结果随后存储在base$ Normalized.PL05数据框中：

> base$Normalized.PL05 <- Vectorize(ChangeFormat)(StormDamageData$Normalized.PL05)

调用包装函数，ChangeFormat：将StormDamageData数据框中的Normalized.CL05作为输入传递。结果随后存储在base$ Normalized.CL05数据框中：

> base$Normalized.CL05 <- Vectorize(ChangeFormat)(StormDamageData$Normalized.CL05)

打印base数据框：head()函数返回基础数据框的前部分。将base数据框作为输入参数传递：

> head(base)

结果如下：

绘制 207 个飓风的标准化成本图：plot()是一个通用函数。base$Normalized.PL05/1e9代表图的x坐标。type="h"代表直方图表示风格。ylim=c(0,155)将 y 轴表示的上限设置为 0（下限）和 155（上限）。x 轴代表损失索引：

> plot(base$Normalized.PL05/1e9, type="h", ylim=c(0,155), main = "207 Hurricanes, Normalized Costs: 1900 - 2005", xlab = "Index of Loss", ylab = "Normalized Costs", col = "red")

结果如下：

步骤 4 - 计算飓风数量

提取每年飓风计数及其频率：基础数据框包含前文所示详情。table()使用base$Year构建每年飓风计数的列联表。结果存储在TestBase数据框中：

> TestBase <- table(base$Year)

打印TestBase数据框的内容：

> TestBase

结果如下：

从 TestBase 数据框中提取年份：names() 函数提取每个年份的名称。as.numeric() 将提取的每个年份名称转换为数值。结果存储在年份数据框中：

> years <- as.numeric(names(TestBase))

打印 years 数据框的内容：

> years

结果如下：

从 TestBase 数据框中提取每年飓风计数的频率：names() 提取每年飓风计数的频率。as.numeric() 将提取的每个飓风计数频率转换为数值。结果存储在频率数据框中：

> frequency <- as.numeric(TestBase)

打印 frequency 数据框的内容：

> frequency

结果如下：

从 TestBase 数据框中提取没有飓风发生年份的飓风计数频率：结果存储在 years0frequency 数据框中：

> years0frequency <- (1900:2005)[which(!(1900:2005)%in%years)]

打印 years0frequency 数据框的内容：

> years0frequency

结果如下：

提取每年飓风的全部计数。结果存储在 StormDamageData 数据框中：

> StormDamageData <- data.frame(years=c(years, years0frequency), frequency=c(frequency, rep(0,length(years0frequency))))

打印 StormDamageData 数据框。head() 函数返回 StormDamageData 数据框的前部分。将 StormDamageData 数据框作为输入参数传递：

> head(StormDamageData) 

结果如下：

绘制 1900 年到 2005 年每年飓风的年份和频率计数：plot() 是一个通用函数。years 代表图表的 x 坐标，而 frequency 代表图表的 y 坐标。type="h" 代表直方图表示风格：

> plot(years, frequency, type="h", main = "Frequency of Hurricanes: 1900 - 2005", xlab = "Time (Years)", ylab = "Annual Frequency", col = "red")

结果如下：

计算从 1900 年到 2005 年所有年份的飓风数量平均值：

> mean(StormDamageData$frequency)

结果如下：

平均每年有两次飓风。

步骤 5 - 构建预测模型

让我们找出飓风发生频率中可能存在的线性趋势。使用 glm() 函数拟合广义线性模型。frequency~years 定义了公式。data = StormDamageData 定义了公式的数据集。family=poisson(link="identity") 函数表示泊松分布。

使用 lm() 函数拟合线性模型。frequency~years 定义了公式。data = StormDamageData 定义了公式的数据集。使用以下命令：

> LinearTrend <- glm(frequency~years, data = StormDamageData, family=poisson(link="identity"), start=lm(frequency~years, data = StormDamageData)$coefficients)

打印 LinearTrend 的详细信息：

> LinearTrend

结果如下：

探索飓风发生频率中可能存在的指数趋势：

使用 glm() 函数拟合广义线性模型。frequency~years 定义公式。data = StormDamageData 定义用于公式的数据集。family=poisson(link="identity") 函数表示泊松分布。我们通过以下命令实现：

> ExpTrend <- glm(frequency~years, data=StormDamageData, family = poisson(link="log"))

打印 ExpTrend 的详细信息：

> ExpTrend

结果如下：

绘制 1900 年至 2005 年间每年飓风的年份和频率计数：plot() 是一个通用函数。年份代表图形的 x 坐标，而频率代表图形的 y 坐标。type="h" 代表直方图表示风格。ylim=c(0,6) 函数将 y 轴表示的上限设置为 0（下限）和 6（上限）：

> plot(years, frequency, type='h', ylim=c(0,6), main = "No. of Major Hurricanes Predicted for 2014", xlim=c(1900,2020))

结果如下：

根据 2014 年的指数趋势预测趋势：使用 predict() 函数根据线性模型对象预测值。ExpTrend 代表从 lm 类继承的对象。newdata = data.frame(years=1890:2030) 函数代表用于查找预测变量的数据框：

> cpred1 <- predict(ExpTrend, newdata = data.frame(years=1890:2030), type="response")

打印 cpred1 的详细信息：

> cpred1

结果如下：

将 cpred1 的点与线段连接：lines() 是一个通用函数，它将 cpred1 数据框的值作为 y 轴的坐标，并用线段连接相应的点。1890:2030 代表 x 轴：

> lines(1890:2030,cpred1,col="blue")

结果如下：

根据 2014 年的线性趋势预测趋势：使用 predict() 函数根据线性模型对象预测值。LinearTrend 代表从 lm 类继承的对象。newdata = data.frame(years=1890:2030) 函数代表用于查找预测变量的数据框：

> cpred0 <- predict(LinearTrend, newdata=data.frame(years=1890:2030), type="response")

打印 cpred0 的详细信息：

> cpred0

结果如下：

将 cpred0 的点与线段连接：lines() 是一个通用函数，它将 cpred0 数据框的值作为 y- 轴的坐标，并用线段连接相应的点。1890:2030 代表 x-轴：

> lines(1890:2030, cpred0, col="red"))

结果如下：

绘制平均值：abline() 作为函数使用平均值绘制直线，该平均值是 StormDamageData$frequency 的 1.95283。h = mean(StormDamageData$frequency) 是水平线的 y 值：

> abline(h = mean(StormDamageData$frequency), col="black")

结果如下：

将 StormDamageData$frequency、cpred0 和 cpred1 的数据框值合并为平均值：

> predictions <- cbind(constant = mean(StormDamageData$frequency), linear = cpred0[126], exponential=cpred1[126])

打印预测的详细信息：

> predictions

结果如下：

在 2014 年的图表上绘制预测点的位置：

> points(rep((1890:2030)[126],3), predictions, col=c("black","red","blue"), pch=19)

结果如下：

重要的是要注意，通过改变预测模型，保费将会发生变化。在平坦的预测中，不到两个（主要）飓风，但以指数趋势，则超过四个（主要）飓风。

步骤 6 - 计算再保险合同的纯保费

现在我们找到一个合适的模型来计算具有免赔额和有限保额的再保险合同的保费。使用希尔尾指数估计器估计尾指数如下。hill() 是用于估计重尾数据的指数尾部的函数，base$Normalized.PL05：

> hill(base$Normalized.PL05)

结果如下：

前面的图示显示，主要飓风的成本具有重尾分布。

将帕累托模型的损失阈值设置为 5 亿以上如下：

> threshold <- .5

我们使用以下命令返回一个表示广义帕累托模型对超过阈值的拟合的 gpd 类对象。数据集由 base$Normalized.PL05/1e9/20 表示：

 > gpd.PL <- gpd(base$Normalized.PL05/1e9/20, threshold)$par.ests

结果如下：

计算超过阈值 0.5 的 base$Normalized.CL05/1e9/20 数据框的均值：

> mean(base$Normalized.CL05/1e9/20> .5)

结果如下：

假设损失超过 5 亿，我们现在可以计算再保险合同的期望值：

 > ExpectedValue <- function(yinf,ysup,xi,beta){
 + as.numeric(integrate(function(x) (x-yinf) * dgpd(x,xi,mu=threshold,beta),
 + lower=yinf,upper=ysup)$value +
 + (1-pgpd(ysup,xi,mu=threshold,beta))*(ysup-yinf))
 + }

按如下方式找到预测数据框的均值：

> predictions[1]

结果如下：

计算超过阈值 0.5 的 base$Normalized.PL05/1e9/20 数据框的均值：

> mean(base$Normalized.PL05/1e9/20>.5)

结果如下：

这表明每个飓风有 12.5% 的可能性会使保险公司损失超过 5 亿。

计算再保险合同的期望值：

> ExpectedValue(2,6,gpd.PL[1],gpd.PL[2])*1e3

结果如下：

这表明再保险公司预期的还款金额约为 3.309865 亿。

计算再保险合同的保费：

> predictions[1] * mean(base$Normalized.PL05/1e9/20> .5) * ExpectedValue(2, 6, gpd.PL[1], gpd.PL[2]) * 1e3

结果如下：

第十四章. 案例研究 - 电力消费预测

简介

电力是唯一一种生产和消费同时进行的商品；因此，在电力市场必须始终保持供应和消费之间的完美平衡。对于任何国家来说，预测电力消费都是国家利益所在，因为电力是能源的关键来源。可靠的能源消费、生产和分配预测符合稳定和长期的政策。规模经济、关注环境问题、监管要求、良好的公众形象，以及通货膨胀、能源价格快速上涨、替代燃料和技术的出现、生活方式的变化等，都产生了使用建模技术的需求，这些技术可以捕捉价格、收入、人口、技术以及其他经济、人口、政策和技术变量的影响。

低估可能导致产能利用率不足，这会导致服务质量下降，包括局部停电，甚至停电。而另一方面，高估可能导致授权一个可能在未来几年内不需要的工厂。要求是确保投资的最佳时机，这是一个长期考虑，合理化定价结构并设计需求管理计划，以满足短期或中期需求的特点。预测进一步推动各种投资、建设和保护计划。

准备工作

为了进行电力消费预测，我们将使用一个收集于智能电表数据的数据集，该数据集按四个位于不同行业的行业进行时间序列汇总。

第 1 步 - 收集和描述数据

应使用名为 DT_4_ind 的数据集。数值变量如下：

• value

非数值变量如下：

• date_time

• week

• date

• type

如何操作...

让我们深入了解。

第 2 步 - 探索数据

以下包需要在第一步加载：

> install.packages("feather")
> install.packages("data.table")
> install.packages("ggplot2")
> install.packages("plotly")
> install.packages("animation")
> library(feather)
> library(data.table)
> library(ggplot2)
> library(plotly)
> library(animation)

注意

版本信息：本页面的代码在 R 版本 3.2.2 中进行了测试

让我们探索数据并了解变量之间的关系。

检查对象是否为 as.data.table()：数据框的二进制列式序列化使用 feather 实现。为了方便在不同数据分析语言之间共享、读取和写入数据，使用 feather。使用 read_feather() 函数读取 feather 文件。

我们将首先导入 DT_4_ind 数据集。我们将把数据保存到 AggData 数据框中：

> AggData <- as.data.table(read_feather("d:/DT_4_ind"))

探索 AggData 数据框的内部结构：str() 函数显示数据框的内部结构。AggData 作为 R 对象传递给 str() 函数：

> str(AggData)

结果如下：

打印 AggData 数据框。head() 函数返回基本数据框的前部分。将 AggData 数据框作为输入参数传递：

> head(AggData)

结果如下：

绘制按行业汇总的电力消耗时间序列数据。

ggplot() 函数声明用于图形的数据框，并指定在整个图形中要共同使用的绘图美学集。data = AggData 是用于绘图的数据库集，而 aes() 描述了数据中的变量如何映射到视觉属性。geom_line() 生成尝试连接所有观察值的单条线：

    > ggplot(data = AggData, aes(x = date, y = value)) +
+ geom_line() + 
    + facet_grid(type ~ ., scales = "free_y") +
    + theme(panel.border = element_blank(),
    + panel.background = element_blank(),
    + panel.grid.minor = element_line(colour = "grey90"),
    + panel.grid.major = element_line(colour = "green"),
    + panel.grid.major.x = element_line(colour = "red"),
    + axis.text = element_text(size = 10),
    + axis.title = element_text(size = 12, face = "bold"),
    + strip.text = element_text(size = 9, face = "bold")) +
    + labs(title = "Electricity Consumption - Industry", x = "Date", y = "Load (kW)")

结果如下：

注意

重要的一点是，与其它行业相比，食品销售与储存行业的消费在假日期间变化不大。

步骤 3 - 时间序列 - 回归分析

回归模型如下所示：

变量（输入）有两种类型的季节性虚拟变量--每日 和每周 。 是时间 i 时的电力消耗，其中 是要估计的回归系数。

打印 AggData 数据框的内容：

> AggData

结果如下：

将工作日的字符转换为整数：使用 as.factor() 函数将向量编码为因子。as.integer() 函数创建 AggData[, week] 的整数类型对象：

> AggData[, week_num := as.integer(as.factor(AggData[, week]))]

打印更改后的 AggData 数据框内容：

> AggData

结果如下：

使用以下方法从 AggData 数据框中提取唯一的行业类型：

 > n_type <- unique(AggData[, type]) 

打印更改后的数据框 n_type 内容：

 > n_type 

结果如下：

使用以下方法从 AggData 数据框中提取唯一日期：

 > n_date <- unique(AggData[, date]) 

使用以下方法从 AggData 数据框中提取唯一的工作日：

 > n_weekdays <- unique(AggData[, week]) 

使用以下方法设置 period 值：

 > period <- 48 

在样本数据集上执行回归分析。

我们在两周的时间内提取教育（学校）建筑。结果存储在 data_reg 数据框中。n_type[2] 代表教育建筑，而 n_date[57:70] 表示两周的时间段：

 > data_reg <- AggData[(type == n_type[2] & date %in% n_date[57:70])] 

打印更改后的 data_reg 数据框内容：

 > data_reg 

结果如下：

在 2 周期间（2 月 27 日至 3 月 12 日）绘制教育（学校建筑）样本数据集：

ggplot() 函数声明了图形的输入数据框并指定了在整个图形中要通用的绘图美学集。data_reg 是用于绘图的数据库，而 aes() 描述了数据中的变量如何映射到视觉属性。geom_line() 生成单条线，试图连接所有观测值：

    > ggplot(data_reg, aes(date_time, value)) +
    + geom_line() +
    + theme(panel.border = element_blank(),
    + panel.background = element_blank(),
    + panel.grid.minor = element_line(colour = "grey90"),
    + panel.grid.major = element_line(colour = "green"),
    + panel.grid.major.x = element_line(colour = "red"),
    + axis.text = element_text(size = 10),
+ axis.title = element_text(size = 12, face = "bold")) 
    + labs(title = "Regression Analysis - Education Buildings", x = "Date", y = "Load (kW)")

结果如下：

从 data_reg 数据框中提取行数：

 > N <- nrow(data_reg) 

计算训练集中的天数：

 > trainset_window <- N / period 

创建独立的季节性虚拟变量--每日 和每周 。每日季节性值从 1,.....period, 1,.......period 中提取 48 个每日变量的向量。每周值从 week_num 中提取。然后将结果存储在一个向量 matrix_train 中：

 > matrix_train <- data.table(Load = data_reg[, value], Daily = as.factor(rep(1:period, trainset_window)), Weekly = as.factor(data_reg[, week_num])) 

在更改后打印 matrix_train 数据框的内容：

 > matrix_train 

结果如下：

创建线性模型。lm() 函数拟合线性模型：Load ~ 0 + . 是公式。由于 lm() 自动添加到线性模型的截距，我们将其现在定义为 0。data = matrix_train 定义了包含数据的数据框：

 > linear_model_1 <- lm(Load ~ 0 + ., data = matrix_train) 

在更改后打印 linear_model_1 数据框的内容：

 > linear_model_1 

结果如下：

生成模型 linear_model_1 的结果摘要：

> summary_1 <- summary(linear_model_1)

在更改后打印 summary_1 数据框的内容：

 > summary_1 

结果如下：

使用 summary_1 数据框中的 r.squared 属性提取决定系数：

> paste("R-squared: ", round(summary_1$r.squared, 3), ", p-value of F test: ", 1-pf(summary_1$fstatistic[1], summary_1$fstatistic[2], summary_1$fstatistic[3]))

从 data_reg 列表创建一个 data.table：

 > datas <- rbindlist(list(data_reg[, .(value, date_time)], data.table(value = linear_model_1$fitted.values, data_time = data_reg[, date_time]))) 

在更改后打印 datas 数据框的内容：

 > datas 

结果如下：

绘制 linear_model_1 的拟合值。

data = datas 是用于绘图的数据库，而 aes() 描述了数据中的变量如何映射到视觉属性。geom_line() 生成单条线，试图连接所有观测值：

 > ggplot(data = datas, aes(date_time, value, group = type, colour = type)) + geom_line(size = 0.8) + theme_bw() + 
 + labs(x = "Time", y = "Load (kW)", title = "Fit from Multiple Linear Regression") 

结果如下：

绘制拟合值与残差值的关系图。

data 是用于绘图的数据库，而 aes() 描述了数据中的变量如何映射到视觉属性：

> ggplot(data = data.table(Fitted_values =
linear_model_2$fitted.values, Residuals = linear_model_2$residuals),
aes(Fitted_values, Residuals)) + geom_point(size = 1.7)
+ geom_hline(yintercept = 0, color = "red", size = 1) +
+ labs(title = "Fitted values vs Residuals")

结果如下：

函数首先给出线性模型的标准化残差。然后计算1Q和4Q线。接着，从正态分布生成分位数分布。然后计算斜率和截距，并将其绘制出来：

    > ggQQ <- function(lm) {
    # extracting residuals from the fit
    + d <- data.frame(std.resid = rstandard(lm))
    # calculate 1Q, 4Q line
    + y <- quantile(d$std.resid[!is.na(d$std.resid)], c(0.25, 0.75))
    # calculate 1Q, 4Q line
    + x <- qnorm(c(0.25, 0.75))
    + slope <- diff(y)/diff(x)
    + int <- y[1L] - slope * x[1L]
+ 
    + p <- ggplot(data = d, aes(sample = std.resid)) +
+ stat_qq(shape = 1, size = 3) + 
+ labs(title = "Normal Q-Q", 
+ x = "Theoretical Quantiles", 
+ y = "Standardized Residuals") + 
    + geom_abline(slope = slope, intercept = int, linetype = "dashed",
+ size = 1, col = "firebrick1") 
    + return(p)
    + }

我们可以使用以下命令绘制 Q-Q 图：

 > ggQQ(linear_model_1) 

结果如下：

如清晰可见，点不在红色线上，它们不正常。由于周变量的估计系数，白天的测量值不断移动，但白天的行为没有被捕捉到。我们需要捕捉这种行为，因为周末的行为尤其不同。

步骤 4 - 时间序列 - 改进回归分析

创建线性模型：lm()函数拟合线性模型。Load ~ 0 + Daily + Weekly + Daily:Weekly是新的公式。由于lm()自动添加到线性模型的截距，我们将其现在定义为0。data = matrix_train定义了包含数据的数据框：

> linear_model_2 <- lm(Load ~ 0 + Daily + Weekly + Daily:Weekly, data = matrix_train)

在更改后打印linear_model_2数据框的内容：

 > linear_model_2 

结果如下：

比较来自linear_model_1和linear_model_2模型摘要的 R-squared 值：

> c(Previous = summary(linear_model_1)$r.squared, New = summary(linear_model_2)$r.squared)

结果如下：

第二个模型的 R-squared 值有显著提高。

图形比较linear_model_1和linear_model_2模型的残差。

 > ggplot(data.table(Residuals = c(linear_model_1$residuals, linear_model_2$residuals), Type = c(rep("Multiple Linear Reg - simple", nrow(data_reg)), rep("Multiple Linear Reg with interactions", nrow(data_reg)))), aes(Type, Residuals, fill = Type)) + geom_boxplot()
 > ggplotly()

结果如下：

linear_model_1的残差细节。

结果如下：

linear_model_2的残差细节。

结果如下：

从data_reg和linear_model_2的列表中创建一个data.table：

 > datas <- rbindlist(list(data_reg[, .(value, date_time)], data.table(value = linear_model_2$fitted.values, data_time = data_reg[, date_time]))) 

在更改后打印datas数据框的内容：

 > datas 

结果如下：

向datas添加Real和Fitted列：

 > datas[, type := rep(c("Real", "Fitted"), each = nrow(data_reg))] 

在更改后打印datas数据框的内容：

 > datas 

结果如下：

绘制linear_model_2的拟合值。

data = datas是用于绘图的数据库集，而aes()描述了数据中的变量如何映射到视觉属性。geom_line()生成试图连接所有观察值的单一线条：

 > ggplot(data = datas, aes(date_time, value, group = type, colour =
type)) + geom_line(size = 0.8) + theme_bw() +
+ labs(x = "Time", y = "Load (kW)", title = "Fit from Multiple Linear
Regression")

结果如下：

与之前 linear_model_1 的绘图相比，拟合值和实际值非常接近。

绘制拟合值与残差值的关系图。Data 是用于绘图的数据库，而 aes() 描述了数据中的变量如何映射到视觉属性：

 > ggplot(data = data.table(Fitted_values = linear_model_2$fitted.values, Residuals = linear_model_2$residuals), aes(Fitted_values, Residuals)) + geom_point(size = 1.7) 
 + geom_hline(yintercept = 0, color = "red", size = 1) + 
 + labs(title = "Fitted values vs Residuals") 

结果如下：

与之前 linear_model_1 的绘图相比，这些图看起来更接近残差线。

我们可以使用以下方式绘制 Q-Q 图：

 > ggQQ(linear_model_2) 

结果如下：

步骤 5 - 构建预测模型

我们可以定义一个函数来返回一周前预测的预测结果。输入参数是 data 和 set_of_date：

    > predWeekReg <- function(data, set_of_date){
    + #creating the dataset by dates
+ data_train <- data[date %in% set_of_date] 
    + N <- nrow(data_train)
    +
    + # number of days in the train set
    + window <- N / period # number of days in the train set
    +
    + #1, ..., period, 1, ..., period - daily season periods
    + #feature "week_num"- weekly season
    + matrix_train <- data.table(Load = data_train[, value],
    + Daily = as.factor(rep(1:period, window)),
    + Weekly = as.factor(data_train[, week_num]))
    +
    + #creating linear model.
    + # formula - Load ~ 0 + Daily + Weekly + Daily:Weekly
    + # dataset - data = matrix_train
    + lm_m <- lm(Load ~ 0 + Daily + Weekly + Daily:Weekly, data = matrix_train)
+ 
    + #forecast of one week ahead
    + pred_week <- predict(lm_m, matrix_train[1:(7*period), -1, with = FALSE])
    + return(as.vector(pred_week))
    + }

定义评估预测的平均绝对百分比误差：

 > mape <- function(real, pred){
 + return(100 * mean(abs((real - pred)/real)))
 + }

将训练集长度设置为 2 周，因此减去 2。将生成 50 周的预测。使用滑动窗口方法进行训练预测，为每种行业进行预测：

> n_weeks <- floor(length(n_date)/7) - 2

打印周数：

> n_weeks

结果如下：

计算每种行业一周前预测的预测结果。

调用函数返回 AggData 商业地产和数据集一周前预测的预测结果：

 > lm_pred_weeks_1 <- sapply(0:(n_weeks-1), function(i)
 + predWeekReg(AggData[type == n_type[1]], n_date[((i*7)+1):((i*7)+7*2)]))

调用函数返回 AggData - 教育和日期集一周前预测的预测结果：

 > lm_pred_weeks_2 <- sapply(0:(n_weeks-1), function(i)
 + predWeekReg(AggData[type == n_type[2]], n_date[((i*7)+1):((i*7)+7*2)]))

调用函数返回 AggData 食品和销售以及日期集一周前预测的预测结果：

 > lm_pred_weeks_3 <- sapply(0:(n_weeks-1), function(i)
 + predWeekReg(AggData[type == n_type[3]], n_date[((i*7)+1):((i*7)+7*2)]))

调用函数返回 AggData 照明行业和日期集一周前预测的预测结果：

 > lm_pred_weeks_4 <- sapply(0:(n_weeks-1), function(i)
 + predWeekReg(AggData[type == n_type[4]], n_date[((i*7)+1):((i*7)+7*2)]))

计算每种行业的平均绝对百分比误差以评估预测。调用函数返回平均绝对百分比。计算评估 AggData 照明行业和日期集预测的误差：

 > lm_err_mape_1 <- sapply(0:(n_weeks-1), function(i)
 + mape(AggData[(type == n_type[1] & date %in% n_date[(15+(i*7)):(21+(i*7))]), value],
 + lm_pred_weeks_1[, i+1]))

打印 lm_err_mape_1 数据框：

> lm_err_mape_1

结果如下：

调用函数返回评估 AggData 教育和日期集预测的平均绝对百分比误差：

 > lm_err_mape_2 <- sapply(0:(n_weeks-1), function(i)
 + mape(AggData[(type == n_type[2] & date %in% n_date[(15+(i*7)):(21+(i*7))]), value],
 + lm_pred_weeks_2[, i+1]))

打印 lm_err_mape_2 数据框：

> lm_err_mape_2

结果如下：

调用函数返回评估 AggData 食品和销售以及日期集预测的平均绝对百分比误差：

 > lm_err_mape_3 <- sapply(0:(n_weeks-1), function(i)
 + mape(AggData[(type == n_type[3] & date %in% n_date[(15+(i*7)):(21+(i*7))]), value],
 + lm_pred_weeks_3[, i+1]))

打印 lm_err_mape_3 数据框：

> lm_err_mape_3

结果如下：

调用函数返回评估 AggData 照明行业和日期集预测的平均绝对百分比误差：

 > lm_err_mape_4 <- sapply(0:(n_weeks-1), function(i)
 + mape(AggData[(type == n_type[4] & date %in% n_date[(15+(i*7)):(21+(i*7))]), value],
 + lm_pred_weeks_4[, i+1]))

打印 lm_err_mape_4data 数据框：

> lm_err_mape_4

结果如下：

步骤 6 - 绘制一年的预测图

绘制结果：

注意

您需要安装 ImageMagick-7.0.4-Q16 以使 saveGIF 功能正常工作。

    > datas <- data.table(value = c(as.vector(lm_pred_weeks_1),
 AggData[(type == n_type[1]) & (date %in% n_date[-c(1:14,365)]), value]),
    date_time = c(rep(AggData[-c(1:(14*48), (17473:nrow(AggData))), date_time], 2)),
    type = c(rep("MLR", nrow(lm_pred_weeks_1)*ncol(lm_pred_weeks_1)),
    rep("Real", nrow(lm_pred_weeks_1)*ncol(lm_pred_weeks_1))),
    week = c(rep(1:50, each = 336), rep(1:50, each = 336)))

    > saveGIF({
    oopt = ani.options(interval = 0.9, nmax = 50)
    for(i in 1:ani.options("nmax")){
    print(ggplot(data = datas[week == i], aes(date_time, value, group = type, colour = type)) +
    geom_line(size = 0.8) +
scale_y_continuous(limits = c(min(datas[, value]), max(datas[, value]))) + 
    theme(panel.border = element_blank(), panel.background = element_blank(),
    panel.grid.minor = element_line(colour = "grey90"),
    panel.grid.major = element_line(colour = "grey90"),
    panel.grid.major.x = element_line(colour = "grey90"),
    title = element_text(size = 15),
    axis.text = element_text(size = 10),
    axis.title = element_text(size = 12, face = "bold")) +
    labs(x = "Time", y = "Load (kW)",
    title = paste("Forecast of MLR (", n_type[1], "); ", "week: ", i, "; MAPE: ",
    round(lm_err_mape_1[i], 2), "%", sep = "")))
    ani.pause()
    }
    }, movie.name = "industry_1.gif", ani.height = 450, ani.width = 750)

结果如下：

前面的结果证明，电力消耗模式是基于外部因素，如假日、天气、物业性质等。消耗模式在本质上是非常随机的。

注意

目标是向读者介绍如何应用多重线性回归来预测双季节时间序列。包含独立变量的交互作用以确保模型的有效性是非常有效的。