# Econometrics vs. Machine Learning with Temporal Patterns

A few months ago, I did publish a (long) post entitled ‘some thoughts on economics, mathematics, econometrics, machine learning, etc‘. In that post, I was discussing possible differences between foundations of econometrics, and machine learning. I wanted to get back today on an important point, related to training/sampling datasets, when we have temporal data.

I was discussing this morning, with a student of the Data Science for Actuaries program, an interesting point related to claim frequency models, for insurance ratemaking. Since the goal is to predict claims frequency (to assess the level of the insurance premium), he suggested to use old data to train the model, and more recent one to test it. The problem is that the model did not incorporate any temporal pattern, and we got surprising results.

Consider here a simple dataset,

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

Claims frequency is driven by a Poisson process, with one covariate, X1, and we assume that the intensity decreases (with an exponential rate). Consider here a standard linear regression, without any time effect

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

We can also compute the empirical annualized claims frequency

```> u=seq(0,1,by=.01)
> v=predict(reg,newdata=data.frame(X1=u,E=1))
> p=function(x){
+   B=B[abs(B\$X1-x)<.1,]
+   sum(B\$Y)/sum(B\$E)
+ }
> vp=Vectorize(p)(seq(.05,.95,by=.1))```

and plot the two curves on the same graph,

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

This is what we usually do in econometrics. In machine learning, and more specifically to assess the quality of the model, and for model selection, it is common to split the dataset in two parts. A training sample, and a validation sample. Consider some randomized training/validation samples, then fit a model on the training sample, and finally use it to get a prediction,

```> idx=sample(1:nrow(B),size=nrow(B)*7/8)
> B_a=B[idx,]
> B_t=B[-idx,]
> reg=glm(Y~X1+offset(log(E)),data=B_a,
+ family=poisson)
> u=seq(0,1,by=.01)
> v=predict(reg,newdata=data.frame(X1=u,E=1))
> p=function(x){
+   B=B_a[abs(B_a\$X1-x)<.1,]
+   sum(B\$Y)/sum(B\$E)
+ }
> vp_a=Vectorize(p)(seq(.05,.95,by=.1))
> plot(seq(.05,.95,by=.1),vp_a,col="blue")
> lines(u,exp(v),lty=2)
> p=function(x){
+   B=B_t[abs(B_t\$X1-x)<.1,]
+   sum(B\$Y)/sum(B\$E)
+ }
> vp_t=Vectorize(p)(seq(.05,.95,by=.1))
> lines(seq(.05,.95,by=.1),vp_t,col="red")```

The blue curve is the prediction on the training sample (as we usually do in econometrics), but then the red curve is the prediction on the testing sample. Here, volatility probably comes from the small size of the testing sample (1 observation out of 8).

Now, what if we use the year as a splitting criteria : we fit a model on old years to fit a model, and we test it on recent years,

```> B_a=subset(B,T<2014)
> B_t=subset(B,T>=2014)
> reg=glm(Y~X1+offset(log(E)),data=B_a,family=poisson)
> u=seq(0,1,by=.01)
> v=predict(reg,newdata=data.frame(X1=u,E=1))
> p=function(x){
+   B=B_a[abs(B_a\$X1-x)<.1,]
+   sum(B\$Y)/sum(B\$E)
+ }
> vp_a=Vectorize(p)(seq(.05,.95,by=.1))
> plot(seq(.05,.95,by=.1),vp_a,col="blue")
> lines(u,exp(v),lty=2)
> p=function(x){
+   B=B_t[abs(B_t\$X1-x)<.1,]
+   sum(B\$Y)/sum(B\$E)
+ }
> vp_t=Vectorize(p)(seq(.05,.95,by=.1))
> lines(seq(.05,.95,by=.1),vp_t,col="red")```

Clearly, we miss something here…

We were looking at such a graph this morning, and it took me some time to understand how training and validation samples were designed, and that there was a possible temporal effect (actually, this morning, it was based on a 3 year training sample, and a 1 year validation sample).

Since there is a temporal pattern, let us capture it. As an econometrician, let me use a regression model

```> reg=glm(Y~X1+T+offset(log(E)),data=B,
+ family=poisson)
> C=coefficients(reg)
> u=seq(1999,2016,by=.1)
> v=exp(-(u-2000)/20-3)
> plot(2000:2015,exp(C[1]+C[3]*(2000:2015)))
> lines(u,v,lty=2,col="red")```

(I focus only on the evolution of the temporal variate on that graph).

Here, we use a linear model, but there are usually no reason to assume linearity. So we might consider splines

```> library(splines)
> reg=glm(Y~X1+bs(T)+offset(log(E)),
+ data=B,family=poisson)
> u=seq(1999,2016,by=.1)
> v=exp(-(u-2000)/20-3)
> v2=predict(reg,newdata=data.frame(X1=0,
+ T=2000:2015,E=1))
> plot(2000:2015,exp(v2),type="b")
> lines(u,v,lty=2,col="red")```

But here again, why should we assume that there is an underlying smooth function? There might be some ruptures… So let us consider a regression on factors

```> reg=glm(Y~0+X1+as.factor(T)+offset(log(E)),
+ data=B,family=poisson)
> C=coefficients(reg)
> u=seq(1999,2016,by=.1)
> v=exp(-(u-2000)/20-3)
> plot(2000:2015,exp(C[2:17]),type="b")
> lines(u,v,lty=2,col="red")```

An alternative might be to consider some more general model, like a regression tree

```> library(rpart)
> reg=rpart(Y~X1+T+offset(log(E)),data=B,
+ method="poisson",cp=1e-4)
> p=function(t){
+   B=B[B\$T==t,]
+   B\$E=1
+   mean(predict(reg,newdata=B))
+ }
> y_m=Vectorize(function(t) p(t))(2000:2015)
> u=seq(1999,2016,by=.1)
> v=exp(-(u-2000)/20-3+.5)
> plot(2000:2015,y_m,ylim=c(.02,.085),type="b")
> lines(u,v,lty=2,col="red")```

Here, it seems that something went wrong. I guess it’s coming from the exposure. So consider a simplier model, on the annualized frequency, and with weights that are related to the exposure

```> reg=rpart(Y/E~X1+T,data=B,weights=B\$E,cp=1e-4)
> p=function(t){
+   B=B[B\$T==t,]
+   B\$E=1
+   mean(predict(reg,newdata=B))
+ }
> y_m=Vectorize(function(t) p(t))(2000:2015)
> u=seq(1999,2016,by=.1)
> v=exp(-(u-2000)/20-3+.5)
> plot(2000:2015,y_m,ylim=c(.02,.085),type="b")
> lines(u,v,lty=2,col="red")```

That was for the econometrician perspective. With a machine learning perspective, consider a training sample (here based on old data) and a validation sample (based on more recent ones)

```> B_a=subset(B,T<2014)
> B_t=subset(B,T>=2014)```

If we consider a model, it is easy to get a prediction on recent years, even if the model was designed to model older ones,

```> reg_a=glm(Y~X1+T+offset(log(E)),
+ data=B_a,family=poisson)
> C=coefficients(reg_a)
> u=seq(1999,2016,by=.1)
> v=exp(-(u-2000)/20-3)
> plot(2000:2015,exp(C[1]+C[3]*c(2000:2013,
+ NA,NA)),type="b")
> lines(u,v,lty=2,col="red")
> points(2014:2015,exp(C[1]+C[3]*2014:2015),
+ pch=19,col="blue")```

But if we use years as factors, things are more complicated.

```> reg_a=glm(Y~0+X1+as.factor(T)+offset(log(E)),
+ data=B_a,family=poisson)
> C=coefficients(reg_a)
> RMSE=function(A){
+   L=exp(C[1]*B_t\$X1+ A[1]*(B_t\$T==2014) + A[2]*(B_t\$T==2015))
+   Y_t=L*B_t\$E
+   sum( (Y_t - B_t\$Y )^2)}
> i=optim(c(.4,.4),RMSE)\$par
> plot(2000:2015,c(exp(C[2:15]),NA,NA),)
> u=seq(1999,2016,by=.1)
> v=exp(-(u-2000)/20-3)
> lines(u,v,lty=2,col="red")
> points(2014:2015,exp(i),pch=19,col="blue")```

becase we need to get a prediction on levels that were not in our training sample. Here, we minimize the RMSE to quantify factor levels for recent years. And the output is not that bad.

So yes, it is possible to get a training dataset on older data, and test it on recent years. But one should be careful, and take into account, properly, temporal patterns.

# From a random generator to a sample function

This week-end, I wrote a post since I had some trouble to generate a sample random sample with R, to reproduce one obtained by a co-author, with SAS (generated using Fishman and Moore (1982) used in function RANUNI). I was lucky since another contributor for that book, Christrophe Dutang, got the anwer to the last question I asked: is it possible to reproduce the random generator ? Yes, we can. And it is quite simple, if you use the appropriate library and the appropriate function,

```> library(randtoolbox)
> a <- 397204094
> b <- 0
> m <- 2^(31)-1
> set.generator(name="congruRand", mod=m, mult=a, incr=b, seed=123)
> runif(10)
[1] 0.7503961 0.3209120 0.1783896 0.9060334 0.3571171
[6] 0.2211140 0.7864383 0.3980819 0.1246652 0.1876858```

If you check in the previous post this is exactly what SAS gave us (and that I could not reproduce by myself). But that was only one part of my problems, since the goal was actually to reproduce indices for a training subsample for credit scoring issues.

I have to admit that I had never though about it before: how should we write a sample function? If values can be replaced, that is fine, we just have to split the unit interval correctly. Like

```> set.seed(95)
> (U=runif(10))
[1] 0.15171155 0.57584087 0.05309844 0.07044485 0.48887914 0.15276707
[7] 0.37405684 0.30006096 0.96997126 0.30373498
> set.seed(95)
> (R=sample(0:99,size=10,replace=TRUE))
[1] 15 57  5  7 48 15 37 30 96 30```

Here, we just truncate from the values obtained from the random generator. And that is just fine. But how do we write a code to sample without replacements? I mean, how do you get that :

```> (S=sample(0:99,size=10,replace=FALSE))
[1] 15 57  5  6 46 14 35 27 89 92```

My initial idea was to use the following technique. The first value is easy to get: we just split the unit interval into 100 subdivision (as before since for the first value, replacement or not, we don’t care) and see in which interval the random value is. And we remove that value from our sample. Then, we split the unit intervall into 99 subdivision, and see in which interval the random value is. It is the 10th? Fine, then our second value is the 10th from our sample (the first value has been removed). Then we split the unit interval in 98 subdivision, etc. The code I wrote to produce that algorithm was the following,

```> mysample1=function(N,unif){
+  n=length(N)
+  size=length(unif)
+  V0=N[trunc(unif[1]*n)+1]
+  N=N[-which(N==V0)]
+  V=V0
+  for(i in 2:length(unif)){
+    V0=N[trunc(unif[i]*(n-i+1))+1]
+    N=N[-which(N==V0)]
+    V=c(V,V0)}
+ return(V)}```

Unfortunetely, I could not reproduce the sample obtained with the R function,

```> mysample1(0:99,unif=U)
[1] 15 58  5  7 49 17 39 31 97 32
> S
[1] 15 57  5  6 46 14 35 27 89 92```

Since Christrophe is an expert on random generators, I did ask him, one more time. And he came up with the following code,

```> mysample2=function(N,unif){
+   integerset=1:length(N)
+   result=rep(NA,length(unif))
+   for(i in 1:length(unif)){
+     intchosen=integerset[ceiling(U[i]*(length(N)-i+1))]
+     integerset[intchosen]=length(integerset)
+     integerset=integerset[-length(integerset)]
+     result[i]=intchosen}
+   return(N[result])}```

which works just fine !

```> mysample2(0:99,unif=U)
[1] 15 57  5  6 46 14 35 27 89 92
> S
[1] 15 57  5  6 46 14 35 27 89 92```

So now, not only can we reproduce random numbers obtained with other software, we can also obtain the same samples indices, with or without replacement ! Thanks Christophe !

[May, 15th] Note that this is note the generator used in SAS, actually. In order to reproduce the sample function of SAS, the algorithm is much more simple, by clearly not that efficicient since we generate a random sample of size 100 (if we have 100 observations), and then, we keep the values associated to the indices of the 10 smallest (if we want a sample of size 10). The code could be

```> mysample3=function(N,unif,size){
+ q=sort(unif)[size]
+ return(N[U<=q])}

> library(randtoolbox)
> a <- 397204094
> b <- 0
> m <- 2^(31)-1
> set.generator(name="congruRand", mod=m, mult=a, incr=b, seed=123) #OK

> U=runif(100)
> mysample3(1:100,U,size=10)
[1] 27 37 47 59 60 71 75 82 84 87```

Thanks Jean-Philippe for the idea (which works).

# Margin of error, and comparing proportions in the same sample

Irecently tried to answer a simple question, asked by @adelaigue. Actually, I thought that the answer would be obvious… but it is a little bit more compexe than what I thought. In a recent survey about elections in Brazil, it was mentionned in a French newspapper that “Mme Rousseff, 62 ans, de 46,8% des intentions de vote et José Serra, 68 ans, de 42,7%” (i.e. proportions obtained from the survey). It is also mentioned that “la marge d’erreur du sondage est de 2,2% ” i.e. the margin of error is 2.2%, which means (for the journalist) that there is a “grande probabilité que les 2 candidats soient à égalité” (there is a “large probability” to have equal proportions).
Usually, in sampling theory, we look at the margin of error of a single proportion. The idea is that the variance of $\widehat{p}$, obtained from a sample of size  is

thus, the standard error is

The standard 95% confidence interval, derived from a Gaussian approximation of the binomial distribution is

The largest value is obtained when p is 1/2, and then we have a worst case confidence interval (an upper bound) which is

So with a margin of error  means that . Hence, with a 5% margin of error, it means that n=400. While 2.2% means that n=2000:
> 1/.022^2
[1] 2066.116
Classically, we compare proportions between two samples: surveys at two different dates, surveys in different regions, surveys paid by two different newpapers, etc. But here, we wish to compare proportions within the same sample. This has been consider in an “old” paper published in 1993 in the American Statistician,

It contains nice figures to illustrate the difference between the standard approach,

and the one we would like to study here.

This point is mentioned in the book by Kish, survey sampling (thanks Benoit for the reference),

Let and denote empirical frequencies we have obtained from the sample, based on  observations. Then since

and

we have

Thus, a natural margin of error on the difference between the two proportion is here

which is here 4 points
> n=2000
> p1=46.8/100
> p2=42.7/100
> 1.96*sqrt((p1+p2)-(p1-p2)^2)/sqrt(n)
[1] 0.04142327
Which is exactly the difference we have here ! Hence, the probability of reaching such a value is quite small (2%)
> s=sqrt(p1*(1-p1)/n+p2*(1-p2)/n+2*p1*p2/n)
> (p1-p2)/s
[1] 1.939972
> 1-pnorm(p1-p2,mean=0,sd=sqrt((p1+p2)-(p1-p2)^2)/sqrt(n))
[1] 0.02619152

Actually, we can compare the three margin of errors we have so far,

• the upper bound
• the “average one”

where

• the more accurate one we just obtained,

where .
> p=seq(0,.5,by=.01)
> ic1=rep(1.96/sqrt(4*n),length(p))
> ic2=1.96*sqrt(p*(1-p))/sqrt(n)
> delta=.01
> ic31=1.96*sqrt(2*p-delta^2)/sqrt(n)
> delta=.2
> ic32=1.96*sqrt(2*p-delta^2)/sqrt(n)
> plot(p,ic32,type=”l”,col=”blue”)
> lines(p,ic31,col=”red”)
> lines(p,ic2)
> lines(p,ic1,lty=2)
So on the graph below, the dotted line is the standard upper bound, the plain line in black being a more accurate one when the probability is  (the x-axis). The red line is the true margin of error with a large difference between candidates (20 points) and the blue line with a small difference (1 point).

Remark: an alternative is to consider a chi-square test, comparering two multinomial distributions, with probabilities  and  where is the average proportion, i.e. 44.75%. Then

i.e.  =3.71
> p=(p1+p2)/2
> (x2=n*((p1-p)^2/p+(p2-p)^2/p))
[1] 3.756425
> 1-pchisq(x2,df=1)
[1] 0.05260495
Under the null hypothesis, should have a chi-square distribution, with one degree of freedom (since the average is fixed here). Here the probability to reach that level is around 5% (which can be compared with the 2% we add before).

So finally, I would think that here, stating that there is a “large probability” is not correct…