# Using convolutions (S3) vs distributions (S4)

Usually, to illustrate the difference between S3 and S4 classes in R, I mention glm (from base) and vglm (from VGAM) that provide similar outputs, but one is based on S3 codes, while the second one is based on S4 codes. Another way to illustrate is to manipulate distributions.

Consider the case where we want to sum (independent) random variables. For instance two lognormal distribution. Let us try to compute the median of the sum.

The distribution function of the sum of two independent (positive) random variables is $F_{S_2}(x)=\int_0^x F_{X_1}(x-y)dF_{X_2}(x)$

pSum2 = function(x) integrate(function(y) plnorm(x-y,1,2)*dlnorm(y,2,1),0,x)$value Let us visualize that cumulative distribution function vx=seq(0.1,50,by=.1) vy=Vectorize(pSum2)(vx) plot(vx,vy,type="l",ylim=c(0,1)) abline(h=.5,lty=2) Let us find an upper bound to compute (in a decent time) quantiles pSum2(350)  0.99195 and then use the uniroot function to inverse that function qSum = function(u) uniroot(function(x) Vectorize(pSum2)(x)-u, interval=c(0,350))$root vu=seq(.01,.99,by=.01) vv=Vectorize(qSum)(vu) The median is here

qSum(.5)  14.155

Why not consider the sum of three (independent) distributions ? Its cumulative distribution function can be writen using our previous function $F_{S_3}(x)=\int_0^x F_{S_2}(x-y)dF_{X_3}(x)$

pSum3 = function(x) integrate(function(y) pSum2(x-y)*dlnorm(y,2,2),0,x)$value If we look at some values we good pSum3(4)  0.015624 pSum3(5) Error in integrate(function(y) plnorm(x - y, 1, 2) * dlnorm(y, 2, 1), : maximum number of subdivisions reached So obviously, there are computational issues here. Let us consider the following alternative expression $F_{S_3}(x)=\int_0^x F_{X_3}(x-y)dF_{S_2}(x)$. Of course, it is necessary here to compute the density of the sum of two variables dSum2 = function(x) integrate(function(y) dlnorm(x-y,1,2)*dlnorm(y,2,1),0,x)$value pSum3 = function(x) integrate(function(y) dlnorm(x-y,2,2)*dSum2(y),0,x)$value Again, let us compute some values pSum3(4)  0.0090285 pSum3(5)  0.01186 This one seems to work quite well. But it is just an illusion. pSum3(9) Error in integrate(function(y) dlnorm(x - y, 1, 2) * dlnorm(y, 2, 1), : maximum number of subdivisions reached Clearly, with those S3-type functions, it wlll be complicated to run computations with 3 variables, or more. Let us consider distributions in the S4-type format of the following package library(distr) X1 = Lnorm(mean=1,sd=2) X2 = Lnorm(mean=2,sd=1) S2 = X1+X2 To compute the median, we simply have to use distr::q(S2)(.5)  14.719 We can also visualize it easily plot(q(S2)) which looks (very) close to what we got, manually. But here, it is also possible to work with the sum of 3 (independent) random variables X3 = Lnorm(mean=2,sd=2) S3 = X1+X2+X3 To compute the median, use distr::q(S3)(.5)  33.208 The function is here plot(q(S3)) # Consommation et production électrique Une version temporaire des slides est en ligne. Pour le TP, on continuera à utiliser les données suivante url="http://freakonometrics.free.fr/elec.RData" download.file(url,"elec.RData") load("elec.RData") k=620 n=nrow(electricite) passe=1:k futur=(k+1):n y=electricite$Load[passe] plot(y,type="l") La semaine passée, nous avions commencé de modéliser la série de température (cette dernière ayant un impact important comme déterminant de la consommation électrique)

y=electricite$Temp[passe] plot(y,type="l") Y=ts(y,start=1996,frequency=52) library(forecast) plot(y[(53:k)-52],y[53:k]) abline(lm(y[53:k]~y[(53:k)-52]),col="red") abline(a=0,b=1,lty=2) La série temporelle est relativement autocorrélée, à l’ordre 52 acf(y,lag=120) model1=auto.arima(Y) acf(residuals(model1),120) On va garder ce modèle en mémoire, et regarder par la suite ses prévisions. Tentons ici un SARIMA, à la main model2=arima(Y,order = c(0,0,0), seasonal = list(order = c(1,0,0))) summary(model2) model2=arima(Y,order = c(0,0,0), seasonal = list(order = c(2,0,0))) summary(model2) acf(residuals(model2),120) model2=arima(Y,order = c(2,0,0), seasonal = list(order = c(1,0,0))) model2=arima(Y,order = c(1,0,0), seasonal = list(order = c(1,1,0))) summary(model2) acf(residuals(model2),120) Tentons ensuite avec une racine unité saisonnière Z=diff(Y,52) model3=arima(Z,order = c(0,0,1), seasonal = list(order = c(0,0,1))) summary(model3) model3=arima(Z,order = c(0,0,2), seasonal = list(order = c(0,0,1))) summary(model3) acf(residuals(model3),120) model3=arima(Y,order = c(0,0,2), seasonal = list(order = c(0,1,1))) summary(model3) On peut ensuite tenter un modèle à la Buys-Ballot reg=lm(Temp~0+as.factor(NumWeek), data=electricite[passe,]) cycle=predict(reg,newdata=electricite[passe,]) plot(y,type="l") lines(cycle,col="red") E=residuals(reg) plot(E,type="l") acf(E,120) model4=arima(E,order = c(1,0,0)) summary(model4) horizon=n-k prev4E=predict(model4,n.ahead = horizon)$pred prev4C=predict(reg,newdata=electricite[futur,]) prev4 =prev4E+prev4C plot(y,type="l",xlim=c(0,n)) abline(v=k,col="red") lines(prev4,col="blue") prev1=predict(model1,n.ahead = horizon)$pred plot(y,type="l",xlim=c(0,n)) abline(v=k,col="red") lines(futur,prev1,col="purple") prev2=predict(model2,n.ahead = horizon)$pred plot(y,type="l",xlim=c(0,n)) abline(v=k,col="red") lines(futur,prev2,col="green") prev3=predict(model3,n.ahead = horizon)$pred plot(y,type="l",xlim=c(0,n)) abline(v=k,col="red") lines(futur,prev3,col="orange") MODELES = data.frame(y1=prev1, y2=prev2,y3=prev3,y4=prev4, y =electricite$Temp[futur]) plot(MODELES$y,type="l") lines(1:111,MODELES$y1,col="purple") lines(1:111,MODELES$y2,col="green") lines(1:111,MODELES$y3,col="orange") lines(1:111,MODELES\$y4,col="blue") 