拓端tecdat|R语言计量经济学与有时间序列模式的机器学习预测

 

今天，我们正在与一位精算数据科学专业的学生讨论，有关保险费率制定的与索赔频率模型有关的观点。由于目标是预测理赔频率（以评估保险费水平），因此他建议使用旧数据来训练该模型，并使用最新数据对其进行测试。问题在于该模型没有包含任何时间模式。

这里考虑一个简单的数据集，

1.  > set.seed(1)
2.  > n=50000
3.  > X1=runif(n)
4.  > T=sample(2000:2015,size=n,replace=TRUE)
5.  > L=exp(-3+X1-(T-2000)/20)
6.  > E=rbeta(n,5,1)
7.  > Y=rpois(n,L*E)
8.  > B=data.frame(Y,X1,L,T,E)

频率由泊松过程驱动，具有一个协变量X1，并且我们假设呈指数速率。在此考虑标准线性回归，没有任何时间因素影响

1.  > reg=glm(Y~X1+offset(log(E)),data=B,
2.  + family=poisson)

我们还可以计算年度经验索赔频率

1.  > u=seq(0,1,by=.01)
2.  > v=predict(re
3.   
4.  > vp=Vectorize(p)(seq(.05,.95,by=.1))

并在同一张图上绘制两条曲线，

1.  > plot(seq(.05,.95,by=.1),vp,type="b")
2.  > lines(u,exp(v),lty=2,col="red")

 

这就是我们通常在计量经济学中所做的。在机器学习中，更具体地说，是评估模型的质量以及进行模型选择，通常将数据集分为两部分。训练样本和验证样本。考虑一些随机的训练/验证样本，然后在训练样本上拟合模型，最后使用它来进行预测，

1.  > idx=sample(1:nrow(B
2.   
3.  > reg=glm(Y~X1+offset(log(E)),data=B_a,
4.  + family=poisson)
5.  > u=seq(0,1,by=.01)
6.  > v=predict(reg,new
7.  $E)
8.  + }
9.  > vp_a=Vectorize(p)(seq(.05,.95,by=.1))
10.  > plot(seq(.05,.95,by=.1),vp_a,col="blue")
11.  > lines(u,exp(v),l
12.  X1-x)<.1,]
13.  + sum(B$Y)/sum(B$E)
14.  + }
15.   
16.  )(seq(.05,.95,by=.1))
17.  > lines(seq(.05,.95,by=.1),vp_t,col="red")

 

蓝色曲线是对训练样本的预测（就像在计量经济学中所做的那样），红色曲线是对测试样本的预测。

现在，如果我们使用年份作为划分标准，我们将旧数据拟合为模型，并在最近几年对其进行测试，

1.  > B_a=subset(
2.  T>=2014)
3.  > reg=glm(Y~X1+offset(l
4.   
5.   
6.  + B=B_a[abs(B_a$X1-x)<.1,]
7.  + sum(B$Y)/sum(B$E)
8.  + }
9.  > vp_a=Vectorize(p)(
10.  y=.1),vp_a,col="blue")
11.  > lines(u,exp(v),lty=2)
12.  > p=function(x){
13.  -x)<.1,]
14.  + sum(B$Y)/sum(B$E)
15.  + }
16.  eq(.05,.95,by=.1))
17.  > lines(seq(.05,.95,b=.1),vp_t,

 

显然，结果不理想。

我花了一些时间来了解训练和验证样本的设计方式对结果产生的影响。

我使用回归模型：

1.   
2.  glm(Y~X1+T+offset(log(E)),data=B,
3.  + family=poisson)
4.   
5.  > u=seq(1999,2016,by=
6.   

 

在这里，我们使用线性模型，但是通常没有理由假设线性。所以我们可以考虑样条

1.   
2.   
3.  > reg=glm(Y~X1+bs(T)+offse
4.   
5.  > u=seq(1999,2016,by=.
6.   
7.  > v2=predict(reg,newdata=
8.   
9.  > plot(2000:2015,exp(v2),ty

 

除了假设存在一个基本的平滑函数，我们可以考虑因子的回归

1.   
2.  as.factor(T)+
3.   
4.  + data=B,family=p
5.  g)
6.  > u=seq(1999,2016,by=.1)
7.  > v=exp(-(u-2000)/20
8.  [2:17]),type="b")

 

另一种选择是考虑一些更通用的模型，例如回归树

1.   
2.   
3.  > reg=rpart(Y~X1+T+offset(log(E)),dat
4.   
5.  > p=function(t){
6.  + B=B[B$T==t,]
7.   
8.  + mean(predict(reg,newdata=B))
9.  + }
10.  2000:2015)
11.  > u=seq(1999,2016,by=.1)
12.   
13.  > plot(2000:2015,y_m,ylim=c(

 

在年化频率上考虑与风险敞口相关的权重 

1.  > reg=rpart(Y/E~X1+T,data
2.  weights=B$E,cp=1
3.   
4.  + B=B[B$T==t,]
5.  + B$E=1
6.   
7.  + }
8.  > y_m=Vectorize(function(t) p(t))(
9.   
10.  > v=exp(-(u-2000)/20-3+.5)
11.  > plot(2000:2015,y_m,ylim=c(.02,.08

 

从机器学习的角度来看，考虑训练样本（基于旧数据）和验证样本（基于较新的样本）

1.  > B_a=subset(B,T<2014)
2.  > B_t=subset(B,T>=2014)

如果我们考虑使用广义线性模型，那么也很容易获得近年来的预测，

1.  > reg_a=glm(Y~X1+T+offset(l
2.   
3.  > C=coefficients(reg_a)
4.  > u=seq(1999,20
5.  /20-3)
6.  +C[3]*c(2000:2013,
7.  + NA,NA)),type="b")
8.  )
9.  > points(2014:2015,exp(C

 

但是，如果我们以年份为因子，我们需要对训练样本中没有的水平进行预测，结果更加复杂。

1.  > reg_a=glm(Y~0+X1+as.factor(T)+offse
2.   
3.  > C=coefficients(reg_a)
4.  2014) + A[2]*(B_t$T==2015))
5.  + Y_t=L*B_t$E
6.   
7.  > i=optim(c(.4,.4),RMSE)$par
8.  > plot(2000:2015,c(exp(C[2:15]
9.   
10.  > lines(u,v,lty=2,col="red")
11.  lue")

我们将RMSE量化近年来的预测水平，输出还不错。

获得旧数据的训练数据集，并在最近几年对其进行测试应该谨慎适当地考虑时间模型。

---