Visualizing uncertainty using Jackknife

Once again, I (re)discovered last week at the Rmetrics conference that old tools can be extremely interesting to illustrate complex ideas, like uncertainty in fnancial markets, and stock prices. For instance a 99.5% quantile: we look for the scenario that occur with a probability of 1 out of 200. Are there nice ways to illustrate that quantity ?

Consider the monthly evolution of the SP500 index over the last 22 years,

> library(quantmod)
> getSymbols('^GSPC', from='1990-01-01')
 "GSPC"
> GSPC = adjustOHLC(GSPC,
+ symbol.name='^GSPC')
> MGSPC = to.monthly(GSPC)
> CLOSE = MGSPC\$GSPC.Close
> plot(CLOSE) It is possible to use Jackknife technique to illustrate uncertainty. The idea, in Jackknife, it to remove one of the observations, and to do that for all observations. More formally, from a sample , we define a (sub)sample where observation as been removed, i.e. . Then, we can study all samples when one observation was removed.

Here, in the context of financial time series, over 270 months, we can wonder what might have been the final value of the index if one observation (i.e. one month) had been removed. It is actually the idea of Jackknife,

> R=diff(log(CLOSE)); R=R[-1]
> n=length(R)
> X=rnorm(n,mean(R),sd(R))
> X=R
> MX=t(matrix(X,n,n))
> MX=exp(MX)
> diag(MX)=1
> SMX=MX
> for(k in 2:n){SMX[,k]=SMX[,k-1]*(MX[,k])}

We can plot the different trajectories of the index, when we remove one month,

> init=as.numeric(CLOSE)
> plot(1:n,init*cumprod(exp(X)),type="l",
+ xlab="",ylab="",col="white")
> for(k in 1:n){lines(0:n,init*c(1,SMX[k,]),
+ col="light blue")}
> lines(0:n,init*c(1,cumprod(exp(X))),lwd=2,
+ col="blue")

This can be used to understand sensitivity, or unccertainty, of financial time series, We can then look closer at the final value of the index, over those 270 scenarios, or we also use a Box-Plot, Here we can clearly see the impact: if we remove one good month, the index ends around 1250, while it reaches 1650 if we remove a bad month. The difference is huge. So instead of talking about volatility (which is actually a complex concept), that Jackknife idea of remove observations might be more intuitive, and much easier to get a first understanding of uncertainty. But those ideas of resampling are great. I will post a nice application soon (but first, I will discuss with some colleagues in Lyon).

“Le prix de l’or bat de nouveaux records”… et alors ?

Histoire de changer un peu, un rapide billet d’économie… bien que l’économie monétaire ne soit pas mon domaine de prédilection, loin de là. Mes compétences en la matière (comme sur beaucoup de notions en économie) se limitent à la lecture (mais aussi une intense relecture) d’Obélix et Compagnie je dois avouer que je suis surpris par les commentaires de certains journalistes économiques, qui nous expliquent que “l’or a atteint des niveaux historiques” (ici, même si ça date un peu) ou que “le prix de l’or bat de nouveaux records” (ce matin cette fois, ). On peut aussi lire ici ou . On peut s’en convaincre en jetant un œil aux graphiques ici, , ou encore. Sans vouloir être méchant, si je transpose ca a une discussion que je pourrais avoir avec mes enfants, ça donnerait: “On a mis 150 euros sur un livret A le jour de votre naissance… c’est incroyable, mais tous les 15 jours, on bat des records: vous êtes de plus en plus riche !” En fait, effectivement, le prix de l’or atteint des records… mais qu’est-ce que je peux faire de ce “prix historique” quand on sait qu’à coté, tout augmente (comme dirait ma grand mère)… On pourrait normaliser par rapport a autre chose, par exemple un indice boursier (car on parle d’investissement), ou comme le font plutôt les économistes, avec quelque chose en rapport avec la masse monétaire. Non ?
Faisons les calculs: les séries du prix de l’or se trouvent ici par exemple (avec des séries journalières pour les années les plus récentes), et pour les données permettant d’actualiser, on pourra aller ici. Pour faire simple, je me suis contenté des données mensuelles en ligne ici, et pour la masse monétaire, M3 serait pas mal (je renvoie ici pour les détails), mais la Réserve Fédérale américaine ne le publie plus depuis 2006. J’aurais bien aimé aller ici, mais c’est payant. Finalement je suis allé , en extrapolant M3 sur les années récentes (de 2006 à 2010) en utilisant le taux de croissance observé sur M2. Les données ne sont pas corrigées des variations saisonnières.

>  OR = read.table("https://perso.univ-rennes1.fr/
arthur.charpentier/prix-or-mensuel.csv",
>  tOR=as.vector(t(OR[,2:13]))
>  VOR=tOR[1:511]

car j’avais mis des 0 pour les mois non observés de 2010 (on les supprime tout simplement). Sinon on retrouve la courbe qui indique qu’effectivement, on a atteint des records…

>    TOR=ts(VOR, frequency = 12, start = c(1968, 1))
>    plot(TOR,col="red") So what ? Regardons un peu notre masse monétaire qui permettra de ramener le prix de l’or dans le passé au prix réel de l’argent,

>    M = read.table("https://perso.univ-rennes1.fr/
arthur.charpentier/M1M2M3.txt")
>    M2=as.numeric(as.character(rev(M\$V4)))
>    M3=as.numeric(as.character(rev(M\$V5)))
>    M2b=M2[560:length(M2)]/M2
>    M3[560:length(M3)]=M3*M2b
>    TM3=ts(M3,frequency=12,start=c(1959,1))
>    plot(TM3,col="blue") Bon, j’ai fait une approximation rapide, mais c’est vraiment histoire de voir ce qui se passe si on normalise le prix de l’or. Ceux qui veulent aller plus loin peuvent améliorer le code. Et en corrigeant pour travailler à quantité d’argent équivalent, on obtient des choses un peu différentes (j’aime les euphémismes)

>    TM3b=ts(M3[(12*9+1):length(M3)],
+    frequency = 12, start = c(1968, 1))
>    TORb=TOR/TM3b
>    plot(TORb,col="purple")
>    abline(h=TORb[length(TORb)],lty=2) On a certes une hausse relativement importante depuis un peu moins de 10 ans du prix de l’or, mais le niveau n’est en rien comparable avec la crise qui s’est observé en 1980 (évoqué ici, , ou encore). Fort heureusement, on ne nous tient pas le même discours sur le PIB, nous expliquant qu’il ne cesse de battre, mois après mois, de nouveaux records ! On notera que ce graphique se rapproche d’autres, en ligne, ici ou
Pour ceux qui ne sont pas convaincu que l’or puisse être relié à une grandeur monétaire, on peut voir l’or comme une alternative à un placement en bourse, et comparer à un indice boursier (par exemple le SP500 pour rester aux États Unis),

table.csv?s=%5EGSPC&a=00&b=1&c=1968&d=08
&e=18&f=2010&g=m&ignore=.csv",
>    tSP=rev(SP\$Close)
>    TSP=ts(tSP, frequency = 12, start = c(1968, 1))
>    plot(TSP,col="red") >    TORc=TOR/TSP
>    plot(TORc,col="purple")
>    abline(h=TORc[length(TORc)],lty=2) Bref, quelle que soit la quantité retenue pour normaliser, on se rend compte que les niveaux atteints ces derniers mois sont très loin de ceux observés en 1980, voire dans les années 80. Pas de quoi en faire la une de la presse….

Nikkei’s past experience vs. SP500 (in euros)

Following Michael’s idea (here), I wanted to go further, based on his intuition (and dataset that he kindly sent me, there). If we consider the two series of Nikkei index and SP500 index in euros, we have to following graph, the code is simply the following (the merging function is simply here to avoid problem with different trading days: since we look at the index and not the return, it is the simplest way to deal with it).

> library(RODBC)
> base = odbcConnectExcel(
+ "https://perso.univ-rennes1.fr/arthur.charpentier/spx_nky_eurusd.xls",
+ readOnly = TRUE)
> series1 = sqlQuery(base,query="select * from [Tabelle1\$A2:B8837]") # SPX
> series2 = sqlQuery(base,query="select * from [Tabelle1\$D2:E8631]") # NKY
> series3 = sqlQuery(base,query="select * from [Tabelle1\$G2:H8945]") # EURUSD
> odbcCloseAll()
> series4=merge(series1,series3)
> series4\$SPEUR=series4\$SPX/series4\$EURUSD
> series5=merge(series4,series2)
> x=(as.Date(series5[,1])-as.Date("01/01/0000","%d/%m/%Y"))/365.25
> yl=range(series5[,4])
> xl=c(1975,2010)
> plot(x,series5[,4],axes=FALSE,xlab="",ylab="",type="l",
+ lwd=3,col="red",xlim=xl,ylim=yl)
> axis(1)
> axis(2, col="red")
> par(new=TRUE)
> yl=range(series5[,5])
> plot(x,series5[,5],axes=FALSE,xlab="",ylab="",type="l",
+ lwd=3,col="blue",xlim=xl,ylim=yl)
> axis(4, col="blue")
> mtext("SP500 in Euros", 2, line=2, col="red", cex=1.2)
> mtext("NKY", 4, line=2, col="blue", cex=1.2)

Those two series series seem to have a similar pattern, so an idea can be translate the SP500 on the left, Interesting isn’t it ? Suppose that we want to forecast (or forsee ?) the SP500 in euro for the next 10 years…

People who enjoy charts would have here a nice tool…

Those two series are extremely correlated, with a correlation of 0.9572,

> X1=series5[2501:n,4]
> X2=series5[1:(n-2500),5]
> cor(X1,X2)
 0.9572484

But are the two series cointegrated (see here, here or therefor material on cointegration) ? Well, using standard procedure, we first have to prove that the two series are integrated. First, let us look at the autocorrelograms,  At first sight, we confirm the economic intuition that those indices should be integrated. Standard tests confirm that intuition,

> acf(X2,lag=1000,col="light green")
> acf(X1,lag=1000,col="light green")
> library(tseries)
Augmented Dickey-Fuller Test
data:  X1
Dickey-Fuller = -1.0768, Lag order = 17, p-value = 0.9264
alternative hypothesis: stationary
Augmented Dickey-Fuller Test
data:  X2
Dickey-Fuller = -1.2905, Lag order = 17, p-value = 0.8788
alternative hypothesis: stationary

But if we want to go further, we have to find the cointegration relationship between the two series. From an heuristic point of view, a linear regression should be a good proxy,

> reg=lm(X1~X2)
> plot(residuals(reg)) > acf(residuals(reg),lag=1000,col="light green") Augmented Dickey-Fuller Test
data:  residuals(reg)
Dickey-Fuller = -5.176, Lag order = 17, p-value = 0.01
alternative hypothesis: stationary
Message d'avis :
In adf.test(residuals(reg)) : p-value smaller than printed p-value
> pp.test(residuals(reg))
Phillips-Perron Unit Root Test
data:  residuals(reg)
Dickey-Fuller Z(alpha) = -46.9775, Truncation lag parameter = 11,
p-value = 0.01
alternative hypothesis: stationary
Message d'avis :
In pp.test(residuals(reg)) : p-value smaller than printed p-value

When we look at the autocorrelation function, it looks like we do have a stationary series.
This idea is – more or less – the idea of Engle-Granger two step procedure. But actually, we can not directly use Dickey-Fuller’s test to see if residuals are integrated. This was proved in Phillips and Ouliaris (1990), who also proposed a test (see e.g. here),

> library(tseries); po.test(cbind(X1,X2))
Phillips-Ouliaris Cointegration Test
data:  cbind(X1, X2)
Phillips-Ouliaris demeaned = -53.1766, Truncation lag parameter = 57,
p-value = 0.01
Message d'avis :
In po.test(cbind(X1, X2)) : p-value smaller than printed p-value
Another similar function can be found in R
> library(urca)
> summary(ca.po(cbind(X1,X2)))
########################################
# Phillips and Ouliaris Unit Root Test #
########################################
Test of type Pu
detrending of series none
Call:
lm(formula = z[, 1] ~ z[, -1] - 1)
Value of test-statistic is: 45.2032
Critical values of Pu are:
10pct    5pct    1pct
critical values 20.3933 25.9711 38.3413

Thus, we has to admit that those series are cointegrated.

Based on that idea, it is possible to model the stationary component, and forecast it for the next ten years, based on the assumption that we know the behavior of one time series. Hence, if we add the confidence interval due to the stationary component uncertainty, we have the following graph, Of course, again, only uncertainty related to the stationary process is considered here….