В курсе о методах резервирования заявок я упоминал об использовании регрессии Пуассона, даже если дополнительные платежи не были целыми числами. Например, мы рассмотрели добавочные треугольники
> source("http://perso.univ-rennes1.fr/arthur.charpentier/bases.R")
> INC=PAID
> INC[,2:6]=PAID[,2:6]-PAID[,1:5]
> INC
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 3209 1163 39 17 7 21
[2,] 3367 1292 37 24 10 NA
[3,] 3871 1474 53 22 NA NA
[4,] 4239 1678 103 NA NA NA
[5,] 4929 1865 NA NA NA NA
[6,] 5217 NA NA NA NA NA
На этих платежах естественно использовать регрессию Пуассона для прогнозирования будущих платежей
> Y=as.vector(INC)
> D=rep(1:6,each=6)
> A=rep(2001:2006,6)
> base=data.frame(Y,D,A)
> reg=glm(Y~as.factor(D)+as.factor(A),data=base,family=poisson(link="log"))
> Yp=predict(reg,type="response",newdata=base)
> matrix(Yp,6,6)
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 3155.6 1202.1 49.8 19.1 8.2 21.0
[2,] 3365.6 1282.0 53.1 20.4 8.7 22.3
[3,] 3863.7 1471.8 60.9 23.4 10.0 25.7
[4,] 4310.0 1641.8 68.0 26.1 11.2 28.6
[5,] 4919.8 1874.1 77.6 29.8 12.8 32.7
[6,] 5217.0 1987.3 82.3 31.6 13.5 34.7
и общая сумма резервов будет
> sum(Yp[is.na(Y)==TRUE]) [1] 2426.985
Здесь платежи были в 000 евро. Что если бы они были в ‘000’000 евро?
> a=1000
> INC/a
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 3.209 1.163 0.039 0.017 0.007 0.021
[2,] 3.367 1.292 0.037 0.024 0.010 NA
[3,] 3.871 1.474 0.053 0.022 NA NA
[4,] 4.239 1.678 0.103 NA NA NA
[5,] 4.929 1.865 NA NA NA NA
[6,] 5.217 NA NA NA NA NA
Мы все еще можем запустить регрессию здесь
> reg=glm((Y/a)~as.factor(D)+as.factor(A),data=base,family=poisson(link="log")) > Yp=predict(reg,type="response",newdata=base) > sum(Yp[is.na(Y)==TRUE])*a [1] 2426.985
и прогноз точно такой же. На самом деле, можно поменять валюту и умножить на любую константу, регрессия Пуассона всегда будет возвращать одно и то же предсказание, если мы используем функцию лог-ссылки,
> homogeneity=function(a=1){
+ reg=glm((Y/a)~as.factor(D)+as.factor(A), data=base,family=poisson(link="log"))
+ Yp=predict(reg,type="response",newdata=base)
+ return(sum(Yp[is.na(Y)==TRUE])*a)
+ }
> Vectorize(homogeneity)(10^(seq(-3,5)))
[1] 2426.985 2426.985 2426.985 2426.985 2426.985 2426.985 2426.985 2426.985 2426.985
Хитрость здесь в том, что нам нравится интерпретация Пуассона. Но GLM просто означает, что мы хотим решить условие первого порядка. Можно явным образом решить условие первого порядка, которое было получено без каких-либо условий, так что значения должны быть целыми числами. Для запуска простого кода перехват должен быть связан с последним значением матрицы, а не с первым.
> base$D=relevel(as.factor(base$D),"6")
> base$A=relevel(as.factor(base$A),"2006")
> reg=glm(Y~as.factor(D)+as.factor(A), data=base,family=poisson(link="log"))
> summary(reg)
Call:
glm(formula = Y ~ as.factor(D) + as.factor(A), family = poisson(link = "log"),
data = base)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.3426 -0.4996 0.0000 0.2770 3.9355
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.54723 0.21921 16.182 < 2e-16 ***
as.factor(D)1 5.01244 0.21877 22.912 < 2e-16 ***
as.factor(D)2 4.04731 0.21896 18.484 < 2e-16 ***
as.factor(D)3 0.86391 0.22827 3.785 0.000154 ***
as.factor(D)4 -0.09254 0.25229 -0.367 0.713754
as.factor(D)5 -0.93717 0.32643 -2.871 0.004092 **
as.factor(A)2001 -0.50271 0.02079 -24.179 < 2e-16 ***
as.factor(A)2002 -0.43831 0.02045 -21.433 < 2e-16 ***
as.factor(A)2003 -0.30029 0.01978 -15.184 < 2e-16 ***
as.factor(A)2004 -0.19096 0.01930 -9.895 < 2e-16 ***
as.factor(A)2005 -0.05864 0.01879 -3.121 0.001799 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 46695.269 on 20 degrees of freedom
Residual deviance: 30.214 on 10 degrees of freedom
(15 observations deleted due to missingness)
AIC: 209.52
Первая идея — запустить градиентный спуск следующим образом (отправной точкой будут коэффициенты из линейной регрессии в журнале наблюдений)
> YNA <- Y
> XNA=matrix(0,length(Y),1+5+5)
> XNA[,1]=rep(1,length(Y))
> for(k in 1:5) XNA[(k-1)*6+1:6,k+1]=k
> u=(1:(length(Y))%%6); u[u==0]=6
> for(k in 1:5) XNA[u==k,k+6]=k
> YnoNA=YNA[is.na(YNA)==FALSE]
> XnoNA=XNA[is.na(YNA)==FALSE,]
> beta=lm(log(YnoNA)~0+XnoNA)$coefficients
> for(s in 1:50){
+ Ypred=exp(XnoNA%*%beta)
+ gradient=t(XnoNA)%*%(YnoNA-Ypred)
+ omega=matrix(0,nrow(XnoNA),nrow(XnoNA));diag(omega)=exp(XnoNA%*%beta)
+ hessienne=-t(XnoNA)%*%omega%*%XnoNA
+ beta=beta-solve(hessienne)%*%gradient}
> beta
[,1]
[1,] 3.54723486
[2,] 5.01244294
[3,] 2.02365553
[4,] 0.28796945
[5,] -0.02313601
[6,] -0.18743467
[7,] -0.50271242
[8,] -0.21915742
[9,] -0.10009587
[10,] -0.04774056
[11,] -0.01172840
Мы не слишком далеки от значений, данных Р. На самом деле, это хорошо, если мы сосредоточимся на прогнозах
> matrix(exp(XNA%*%beta),6,6))
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 3155.6 1202.1 49.8 19.1 8.2 21.0
[2,] 3365.6 1282.0 53.1 20.4 8.7 22.3
[3,] 3863.7 1471.8 60.9 23.4 10.0 25.7
[4,] 4310.0 1641.8 68.0 26.1 11.2 28.6
[5,] 4919.8 1874.1 77.6 29.8 12.8 32.7
[6,] 5217.0 1987.3 82.3 31.6 13.5 34.7
которые именно то, что получено выше. И здесь мы ясно видим, что нет такого предположения, как « объясненная переменная должна быть целым числом ». Также можно помнить, что условие первого порядка совпадает с тем, которое мы имели для взвешенной модели наименьших квадратов. Проблема в том, что весовые коэффициенты являются функцией прогноза. Но используя итерационный алгоритм, мы должны сходиться,
> beta=lm(log(YnoNA)~0+XnoNA)$coefficients
> for(i in 1:50){
+ Ypred=exp(XnoNA%*%beta)
+ z=XnoNA%*%beta+(YnoNA-Ypred)/Ypred
+ REG=lm(z~0+XnoNA,weights=Ypred)
+ beta=REG$coefficients
+ }
>
> beta
XnoNA1 XnoNA2 XnoNA3 XnoNA4 XnoNA5 XnoNA6
3.54723486 5.01244294 2.02365553 0.28796945 -0.02313601 -0.18743467
XnoNA7 XnoNA8 XnoNA9 XnoNA10 XnoNA11 -0.50271242 -0.21915742 -0.10009587 -0.04774056 -0.01172840
которые имеют те же значения, что мы получили ранее. Здесь снова, прогноз такой же, как мы получили из этой так называемой регрессии Пуассона,
> matrix(exp(XNA%*%beta),6,6)
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 3155.6 1202.1 49.8 19.1 8.2 20.9
[2,] 3365.6 1282.0 53.1 20.4 8.7 22.3
[3,] 3863.7 1471.8 60.9 23.4 10.0 25.7
[4,] 4310.0 1641.8 68.0 26.1 11.2 28.6
[5,] 4919.8 1874.1 77.6 29.8 12.8 32.7
[6,] 5217.0 1987.3 82.3 31.6 13.5 34.7
Опять же, это работает просто отлично, потому что GLM представляют собой в основном условия для первых двух моментов, а численные вычисления основаны на условии первого порядка, которое имеет меньше ограничений, чем интерпретация в терминах модели Пуассона.