GRS

may 062013
 
tree_map_export_esp_all_show_2010

Mientras las Administraciones Públicas españolas continúan sin ser capaces de proporcionar unas interfaces homogéneas y modernas para que cualquier analista pueda acceder a los datos disponibles, el mundo se mueve en otra dirección.

Últimamente he sido consciente de la aparición de unas cuantas webs de datos que añaden una capa de depuración y visualización a fuentes ya disponibles. Este enriquecimiento hace más fiable la información de partida, pero especialemente, más accesible y cómoda de utilizar. A continuación las presento y animo a los lectores a explorarlas (todas son gratuítas).

Quandl

La web de Quandl (http://www.quandl.com) proporciona acceso a una amplia selección de datos que proceden de diferentes fuentes en distintos ámbitos: mercados, economía, demografía, cotizaciones, y un largo etcétera. Los datos se pueden descargar y visualizar on-line y dispone de un paquete de funciones de R que facilitan la descarga de datos en dicho entorno de análisis.

Knoema

Knoema (http://knoema.com) tiene una filosofía y cobertura de variables similar a Quandl, pero pone el énfasis en la elaboración de visualizaciones sofisticadas que se pueden incrustar en páginas web. El mecanismo de extracción de datos es menos refinado que en Quandl, envían un email automático con los datos solicitados en formato .xls.

The observatory of economic complexity

Partiendo de los datos brutos de comercio internacional, R. Hausmann y sus discípulos han hecho accesible una riquísima fuente de información estadística, homogeneizando y alargando muchas de las series disponibles (http://atlas.media.mit.edu/). Las impactantes visualizaciones de exportaciones e importaciones proporcionan un resumen compacto del intercambio de bienes entre países y los datos se pueden descargar. Si bien el nivel de desagregación es suficiente para muchos usos, se hace necesario recurrir a COMTRADE para obtener información por país, año y producto.

Paquete WDI para R

El Banco Mundial, bajo la denominación World Development Indicators (http://data.worldbank.org) elabora una base de datos con más de 8000 series temporales comparables entre países que abarcan todo el espectro económico y demográfico. Tanto Quandl como Knoema utilizan ampliamente dicha base de datos y se pueden descargar de la web, pero el paquete WDI de V. Arel-Bundock facilita el acceso directo, selección y carga de datos en R.

mar 012013
 
Captura de pantalla 2013-03-01 a la(s) 13.50.50

He tenido el blog un poco abandonado últimamente, pero no he dejado de enredar. Lo último que he estado desarrollando con R es una aplicación web para la identificación y estimación de modelos ARIMA. La idea surgió cuando, en el primer cuatrimestre, volví a enseñar este tipo de modelos a la primera promoción del Grado de Economía de la Especialidad de Análisis Económico (aunque no llegaron a ver la aplicación).

El paquete shiny de R apareció a finales de 2012 y hace sencillo crear pequeñas aplicaciones a las que el usuario accede desde el navegador. Pueden estar alojadas en un servidor externo o en el propio equipo del usuario. Con dicha herramienta, una vez se adquiere cierta soltura, se puede crear aplicaciones vistosas en relativamente poco tiempo.

Ayer presenté la versión actual (28 de febrero de 2013) de la aplicación en el Grupo de Usuarios de R de Madrid. El video de la presentación se puede ver pulsando aquí. La aplicación en vivo aparece a continuación (es trabajo en curso, así que habrá muchos cambios y puede no funcionar). También se puede acceder a ella en http://glimmer.rstudio.com/grserrano/idts.

Funciona mejor en Chrome que en otros navegadores, pero estaré encantado de recibir comentarios con las combinaciones de navegador y sistema operativo donde se aprecie algún problema. En el iPhone se ve estupendamente, así que ya tengo un ARIMAphone:

idts1idts2

 

 

 

 

 

 

 

sep 042012
 
tabla-isp

Los organismos públicos españoles no ponen sencillo a los analistas el trabajo con los datos que publican (que son menos de los que sería deseable en un país con verdadera transparencia). Cada entidad que publica datos lo hace con su propio formato, que en ocasiones es el de algún programa estadístico comercial (¿para cuándo formatos abiertos?). En algunos casos, como el de este ejemplo, la información se publica como una tabla en una página web, de la que hay que extraer los datos de interés.

Continue reading »

sep 032012
 
Espresso

Actualización 5 de septiembre: Me ha llamado la Decana (Begoña García Greciano) para explicarme que los nuevos precios están aprobados para toda la UCM, no es algo que dependa de nuestra cafetería. Gracias, Begoña.


Si es que tenemos cara de tontos o se nos está poniendo. El gobierno primero nos reduce el sueldo, luego lo congela y después nos lo vuelve a reducir (gobiernos diferentes, pero lo mismo me da). Pero claro, eso denota que somos unos pánfilos y la empresa que gestiona la cafetería se ha dado cuenta de ello y aprovecha: caña al mono, que es de goma.
Continue reading »

sep 032012
 
chiste-politica-1

Quería preparar un ejemplo de web scraping (extracción de datos de páginas web) con R y me decidí por utilizar algún indicador de los que publica periódicamente el Centro de Investigaciones Sociológicas (CIS) en los Barómetros mensuales. Elegí los indicadores de la situación política. Ciertamente, no es que anduviera sobrado de esperanza en lo tocante a la confianza que los ciudadanos vienen depositando en los políticos, pero los datos son para llorar, y es que viene de antiguo:

Continue reading »

may 302012
 
Rlogo-beta

Otro sencillo ejemplo, pero relacionado con el anterior, nos permitirá estudiar la potencia del contraste de Jarque-Bera para distintos tamaños muestrales. La potencia es la probabilidad de rechazar la hipótesis nula cuando es falsa, así que vamos a generar 1000 muestras tamaño 100 de una variable aleatoria U(0,1) y calcularemos si la nula de normalidad se rechaza para cada una de ellas; es decir, vamos a contar el número de veces que el estadístico JB es mayor que 5.99 (el valor crítico de una χ22 para un nivel de significación del α = 0.05).

> set.seed(301)  ## Para hacer el ejemplo reproducible
> nrep  <- 1000
> n     <- 100
> Z     <- matrix(runif(n*nrep, 0, 1), nrow=n)
> jbrechaza <- function(x) {
+     return(jarque.bera.test(x)$statistic > qchisq(0.95, 2))
+ }
> rechaza  <- apply(Z, 2, jbrechaza)
> (potencia <- sum(rechaza)/length(rechaza))
 


[1] 0.549

Incluso con muestras relativamente grandes para lo que se maneja a menudo en econometría, el estadístico de Jarque-Bera sólo es capaz de rechazar la nula de normalidad en una de cada dos ocasiones en que la muestra es uniforme. Como cabe esperar, con muestras menores la potencia se reduce apreciablemente, acercándose a cero con n=20. Puede ir cambiando la variable n para trazar la curva de potencia del contraste para una alternativa U(0,1).

may 292012
 
Rlogo-beta

En R resulta especialmente cómodo programar simulaciones para estudiar propiedades estadísticas de contrastes y estimadores. En este post y los dos siguientes se presentan algunos ejemplos.

Vamos a ilustrar el teorema central del límite, que en una versión sencilla podemos formular del siguiente modo. Sean Xi variables aleatorias independientes y no muy diferentes:


W_n = \sum _{i=1}^n X_i \qquad \Rightarrow \qquad \frac{W_n -    \mathrm{E}(W_n)}{\sqrt{\mathrm{V}(W_n)}} \underset{n \to \infty}{\longrightarrow} N(0,1).

Para una aplicación empírica, se generan 1000 muestras, tamaño 2000, de una población log-normal (supuesto muy frecuente en econometría, donde se toman logaritmos para inducir normalidad), que son las Xi, y las almacenamos en una matriz 2000×1000. Si sumamos los primeros 10 elementos de cada columna, obtenemos una muestra tamaño 1000 de W10, que representamos en un histograma, no se parece a una normal (Gráfico 1(a)).

(a) Histograma de W10

Histograma de W500

(b) Histograma de W500

(c) Histograma de W2000


Gráfico 1: Convergencia de la suma de log-normales a la normal.

Después tomamos los primeros 500 elementos de cada muestra de log normales y los sumamos, tendríamos una muestra de W500. El histograma de W500 empieza a mostrar un comportamiento más simétrico (Gráfico 1(b)). Finalmente, cogemos los 2000 elementos de cada muestra y los sumamos. Ahora disponemos de una muestra tamaño 1000 de la variable aleatoria W2000 (Gráfico 1(c)). El script de R para hacer esta simulación es:

> set.seed(301) ## Para hacer el ejemplo reproducible
> n    <- 2000
> nrep <- 1000
> X    <- matrix(exp(rnorm(n*nrep)), nrow=n)
> W10   <- apply(X[1:10, ], 2, sum)
> W500  <- apply(X[1:500, ], 2, sum)
> W2000 <- apply(X[1:2000, ], 2, sum)
> sapply(data.frame(W10, W500, W2000), jarque.bera.test)[c(1,3),]


          W10      W500         W2000
statistic 7842.408 32.48236     1.754108
p.value   0        8.841903e-08 0.4160067

La función sapply() aplica la función que aparece como segundo argumento a las columnas del data.frame en su primer argumento. De todos los resultados que devuelve la función jarque.bera.test(), en este punto sólo nos interesan dos de ellos: el valor del estadístico y el p-valor, que aparecen en la primera y tercera posición.

Utilizando el contraste de Jarque y Bera, se rechaza la nula para W10 y W500, pero no se puede rechazar normalidad para W2000. Para comparar la tasa de convergencia, puede cambiar la generación de una log-normal (exp(rnorm(n*nrep))) por una uniforme: X <- matrix(runif(n*nrep), nrow=n), para comprobar que a partir de 5 ó 10 uniformes, la suma es aproximadamente normal.

may 172012
 
Rlogo-beta

5.5  Cointegración y modelo de correción de error

Se habla de cointegración entre dos series temporales cuando ambas son no estacionarias (requieren una diferencia para serlo) pero una combinación lineal de ellas sí lo es. En este caso, la combinación lineal estacionaria se puede interpretar como una relación a largo plazo e incluirla retardada en el modelo en diferencias como factor de ajuste o término de corrección de error.

Aunque sin justificación teórica, nos podemos preguntar si el consumo de cigarrillos y la publicidad están cointegradas. Para ello estimamos un modelo en logaritmos de las variables y calculamos los residuos, que deberían ser estacionarios. Para contrastar formalmente la estacionariedad se utiliza el contraste de Dickey-Fuller aumentado (adf.test() del paquete tseries), aplicado a los residuos de la regresión en logaritmos de consumo de tabaco sobre publicidad:

> library(tseries)
> mcoint <- dynlm(log(consumo)~log(publicidad),
+                 data=window(tabdz, start=1940))
> summary(mcoint)


Time series regression with "zooreg" data:
Start = 1940, End = 1978

Call:
dynlm(formula = log(consumo) ~ log(publicidad), data = window(tabdz,
    start = 1940))

Residuals:
     Min       1Q   Median       3Q      Max
-0.56682 -0.01087  0.03205  0.07480  0.16525

Coefficients:
                Estimate Std. Error t value Pr(>|t|)
(Intercept)      7.96733    0.04374 182.138  < 2e-16 ***
log(publicidad)  0.17307    0.04195   4.126 0.000201 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1445 on 37 degrees of freedom
Multiple R-squared: 0.3151,	Adjusted R-squared: 0.2966
F-statistic: 17.02 on 1 and 37 DF,  p-value: 0.0002011
 


> adf.test(as.ts(residuals(mcoint)))  ## Requiere ts, no sirve zoo


	Augmented Dickey-Fuller Test

data:  as.ts(residuals(mcoint))
Dickey-Fuller = -3.3081, Lag order = 3, p-value = 0.08536
alternative hypothesis: stationary

No se puede rechazar la presencia de raíz unitaria en los residuos al 5% de significación (al no poder rechazar la nula), aunque sí al 10%. Si se asume la estacionariedad de los residuos de la regresión anterior, se podrían incluir retardados en la estimación del modelo como un término de corrección de error o relación a largo plazo entre las variables. A continuación, añadimos el término de corrección de error (residuos de la estimación en logaritmos) al data.frame tabz con merge() para, después, estimar los modelos tabm1 y tabm2 incluyendo dicha variable retardada:

> tabdz <- merge(tabdz, err.c=residuals(mcoint))
> tab.ecm1 <- dynlm(d(log(consumo))~d(log(publicidad))+L(err.c),
+                  data=window(tabdz, start=1940))
> tab.ecm2 <- dynlm(d(log(consumo))~L(d(log(consumo)))+L(err.c),
+                  data=window(tabdz, start=1940))
> mtable(tab.ecm1, tab.ecm2)


Calls:
tab.ecm1: dynlm(formula = d(log(consumo)) ~ d(log(publicidad)) + L(err.c),
    data = window(tabdz, start = 1940))
tab.ecm2: dynlm(formula = d(log(consumo)) ~ L(d(log(consumo))) + L(err.c),
    data = window(tabdz, start = 1940))

=======================================
                    tab.ecm1  tab.ecm2
---------------------------------------
(Intercept)          0.021***  0.014*
                    (0.006)   (0.006)
d(log(publicidad))  -0.083
                    (0.075)
L(err.c)            -0.217*** -0.219***
                    (0.042)   (0.054)
L(d(log(consumo)))             0.272*
                              (0.125)
---------------------------------------
R-squared              0.545     0.537
adj. R-squared         0.519     0.509
sigma                  0.033     0.031
F                     20.983    19.688
p                      0.000     0.000
Log-likelihood        77.507    77.476
Deviance               0.038     0.033
AIC                 -147.014  -146.951
BIC                 -140.464  -140.508
N                     38        37
=======================================