Series Temporales

Ejercicio 1

El proceso de modelización de series temporales implica una serie de etapas. Ordene cronológicamente (añadiendo un número de orden 1,2,3,…) cada una de las siguientes etapas en el proceso de Análisis series temporales. Discuta brevemente si lo considera preciso

  1. Representación gráfica de la serie (permea todo el proceso)

  2. Análisis de datos anómalos y tratamiento en su caso

  3. Análisis de datos faltantes e interpolación en su caso

  4. Desestacionalización de la serie

  5. Test de Raices Unitarias

  6. Análisis de la Cointegración a Largo Plazo (sujeto a la disponibilidad de otras series)
    • Diagnosis de los Residuos: Test Ljung-Box (para la regresión de cointegración)
    • Diagnosis de los Residuos: Test BDS (para la regresión de cointegración)

  7. Diferenciación de la Serie para conseguir Estacionariedad

  8. Modelización con modelos ARMA/ARMAX

  9. Análisis de la homocedasticidad de los residuos y modelización ARMA/GARCH

  10. Diagnosis de los Residuos: Test Ljung-Box

  11. Diagnosis de los Residuos: Test BDS

Ejercicio 2

A continuación le mostramos el gráfico de la serie de PIB, datos Brutos, Índice de Volumen. ¿Aprecia a simple vista algún tipo de comportamiento en la serie que requiera su consideración?
PIB - Índice de Volumen

PIB - Índice de Volumen

Se trata de una serie en la que podemos distinguir con mucha claridad sus componentes:

  • Muestra claramente una tendencia ascente, si hacemos el promedio en ventanas temporales veremos que aumenta de forma aproximadamente lineal

  • Muestra además una componente estacional también muy marcada, ya que la estructura de picos parece irse repitiendo a lo largo de la cuurva

No podemos descartar la presencia del efecto calendario pero no es algo reconocible a simple vista sobre la serie.

Ejercicio 3

Indique que transformación realiza cada una de estas sentencias de código R aplicadas sobre la serie de frecuencia trimestral llamada y.ts
  1. Esta línea divide la serie entre si misma avanzada 4 periodos en el tiempo (1 año) y resta uno, de tal forma que nos devuelve la tasa de variación interanual, es decir, el porcentaje de variación de la serie con respecto a sí misma en el pasado ($ = $).
  1. Esta línea hace lo mismo que la línea anterior pero con la serie desfasada un periodo (1 trimestre).
  1. Esta línea devuelve los valores de la serie que están a dos desviaciones estándar o menos de la media. La función scores(…, type=z) normaliza los datos mediante: \(x' = \tfrac{x - \bar{x}}{s_x}\).
  1. Esta línea emplea un filtro de Kalman estacional para rellenar los NANs presentes en la serie, y devuelve la serie sin ellos.
  1. Esta línea lo que hace es promediar el valor de la serie en una ventana trimestral que la recorre de principio a fin, dado que la serie del enunciado es trimestral, no apreciaríamos ningún efecto, tendría más sentido emplear to=“annual” o to=“monthly”.

Ejercicio 4

Explique brevemente que es el componente de Semana santa y calendario de una serie temporal. ¿En qué otros componentes se puede descomponer una serie temporal?

Llamamos componente calendario a aquel comportamiento de la serie debido puramente a efectos relacionados con el calendario, con la forma definida adhoc por nosotros para regular el tiempo, por ejemplo, el efecto del número de días laborales de los meses, las festividades, la no coincidencia entre año natural y lectivo, etc. El componente Semana Santa se puede identificar, por lo tanto, como un efecto calendarío, ya que consite en comportamientos irregulares de la serie en torno a estas fechas, en las que posiblemente exista una disminución en la producción y un aumento en el consumo, por ejemplo. A veces es sútil la diferencia con la componente estacional ya que esta componente recoge un efecto sobre la serie que se repite cada cierto periodo de tiempo y el calendario tiene sus propios periodos, la mejor forma de distinguirlas es pensar si el efecto se mantendría si cambiásemos el calendario, Aparte de las componentes comentadas nos faltan dos, por un lado la tendencia que recoge el movimiento promedio de la serie a lo largo del tiempo, y nos permite identificar un comportamiento general, por el otro el error que se trata de la componente que no somos capaces de explicar de la serie, aquella que se comporta de forma de completamente aleatoria y en la que no identificamos ningún patrón.

Suponga que utiliza la función seas() de la librería de R seasonal para descompner la serie temporal trimestral y.ts
¿Con qué comando podría guardar en otra serie que se llamase yy.ts el componente de ciclo y tendencia de dicha descomposición y_sa?

Para extraer el componente estacional podemos usar seasonal(y_sa) y para el componente de tendencia trend(y_sa), aunque esto funciona mucho mejor aplicado sobre la salida de la función decompose, pese a que esta función aplica un modelo de media móvil más simple.

Ejercicio 5

¿Para qué sirve este test?, ¿cuál es la hipótesis nula y que significa?¿cuál sería el orden de integración de la serie y.ts?

La prueba de Dickey-Fuller trata de determinar si existen raíces unitarias de la serie temporal, lo que implicaría que NO es estacionaria. La hipótesis nula (\(H_0\)) es que si existen tales raices (luego que la serie no es estacionaria) y la rechazaremos cuando el p-valor sea inferior a una significancia establecida (pongamos 0.05). De esta forma el orden de integración de la serie sería 2 (\(I(2)\)) ya que el p valor asociado a las segundas diferencias es 0.01 < 0.05.

Ejercicio 6

¿En qué consiste la cointegración de dos series temporales? Explique la estrategia empírica a seguir para contrastar si dos series temporales están cointegradas o no

Dos series temporales \({X, Y}\) decimos que están cointegradas si ambas poseen el mismo orden de integración \((X \sim I(d)\) e \(Y \sim I(d)\)) y además existen los coeficientes \({\alpha, \beta}\) tal que: \(Z = \alpha X + \beta Y \sim I(0)\), de esta forma nos garantizamos que existe relación y que dicha relación no es espuria. El procedimiento empírico puede consitir en proceder a la comprobación de las condiciones descritas anteriormente o se puede seguir el método de Engle-Granger que consiste en comprobar si el orden de integración de los residuos del modelo de regresión de cointegración es menor al de las series.

Ejercicio 7

Explique en qué consisten los modelos de corrección del error. ¿Por qué se dice que son modelos que combinan el largo plazo con el análisis a corto plazo?

Los modelos de corrección del error emplean la desviación a largo plazo existente entre las series temporales cointegradas y su modelo de regresión de cointegración (tendencia) para definir su dinámica a corto plazo. Se definen mediante la formula: \[ \Delta y_t = \phi_0 - \gamma_1 \hat{\pi}_{t-1} + \sum_{j=1}^J \phi_{1j}\ \Delta x_{t-j} + \sum_{k=1}^K \phi_{2k}\ \Delta y_{t-k} + a_t\\ Donde:\ \hat{\pi}_{t-1} = y_{t-1} - \hat{y}_{t-1} = y_{t-1} - \hat{\beta}_0 - \hat{\beta}_1 x_{t-1} = \text{diferencia serie - tendencia último periodo} \]

Ejercicio 8

8.1

¿Qué es el coeficiente de correlación cruzada entre dos series temporales?. A la vista de este ccf entre la serie ICC.ts.Trim (Indice de Comfianza del Consumidor Trimestral) y la serie TPIB.ts (Tasa de Crecimiento del PIB), ¿diría que la máxima correlación entre ambas series es contemporánea?, ¿hay alguna serie que adelante el comportamiento a la otra?
Correlación Cruzada. ICC - PIB

Correlación Cruzada. ICC - PIB

El coeficiente de correlación cruzada entre dos series temporales es: \[ r_{xy}(T) = \frac{\sigma_{xy}(T)}{\sqrt{\sigma_{xx}(0)\ \sigma_{yy}(0)}}\\ \text{con: }\ \sigma_{xy}(T) = \frac{1}{N-1} \sum_{t=1}^N(x_{t-T} - \bar{x})(y_t - \bar{y}) \] Que es la definición del coeficiente de correlación entre dos variables cualesquiera pero una de ellas la retrasamos \(T\) periodos, de forma que el coeficiente va variando según lo haga el valor de \(T\).

Para las variables de la imagen, vemos que el máximo de la función de correlación cruzada se produce en \(T=-2\), es decir, que el PIB tarda seis meses en reaccionar al comportamiento del ICC, o con mayor lógica, que los consumidores anticipan seis meses el comportamiento del PIB.

8.2

¿Qué es la causalidad en sentido de Granger? Observe los siguientes resultados del test de granger de causalidad entre ICC.ts.Trim y la serie TPIB.ts, podría afirmar que es el ICC (Indice de confianza del consumidor la serie que adelanta al Crecimiento del PIB, os es justo al contrario?¿por qué?

La causalidad den sentido Granger alude exclusivamente a la predicción, y se emplea para determinar si una serie es útil de cara a predecir el comportamiento de otra. Esta causalidad no pretende captar la presencia de efectos comunes y es incapaz además de captar relaciones instantáneas o no lineales, no obstante es un concepto fuertemente asentando en estadística dada su simplicidad computacional.

A la vista de los resultados del test, que determina si la serie X (primer argumento) incorpora información a la predicción de la serie Y (segundo argumento), frente al uso exclusivamente del pasado de la serie Y. Mirando los valores de la significancia, podemos decir que claramente el ICC es útil a la hora de describir el comportamiento del PIB, ya que como vimos en el gráfico ccf va por delante 6 meses (2 lags).

Ejercicio 9

¿Qué es un modelo ARIMA?. ¿Podría escribir la forma analítica que tendría un proceso ARIMA (2,1,1)?

Los modelos AutoRegresive Integrated Moving Average (ARIMA) tratan de explicar el comportamiento de una serie temporal mediante su pasado. Los modelos ARIMA para series temporales sin componente estacional suelen denotarse ARIMA(p, d, q) donde la p es el número de lags involucrados, d es el orden de diferenciación requerido por la serie hasta ser estacionaria y q es el orden de la regresión lineal que explica el error.

La forma más elegante de describirlo matemáticamente es partiendo de un modelo ARMA(p, q) (sin parte integrada): \[ X_t = AR(p) + MA(q) = \alpha_1 X_{t-1} + \alpha_2 X_{t-2} + \ldots + \alpha_p X_{t-p} + \epsilon_t + \theta_1 \epsilon_{t-1} + \theta_2 \epsilon_{t-2} + \ldots + \theta_q \epsilon_{t-q} \] Si llamamos \(L^k\) al operador lag que actúa tal que: \(L^k X_t = X_{t-k}\) podemos escribir el modelo ARMA(p, q) como: \[ X_t = \left(\sum_{i=1}^p \alpha_i L^i\right) X_t + \left(\sum_{i=1}^q \theta_i L^i\right) \epsilon_t + \epsilon_t \longrightarrow \left(1 - \sum_{i=1}^p \alpha_i L^i\right) X_t = \left(1 + \sum_{i=1}^q \theta_i L^i\right) \epsilon_t \] Finalmente para incorporar la parte integrada (\(I(d)\)) debemos suponer que el operador \(\left(1 - \sum_{i=1}^p \alpha_i L^i\right)\) tiene como factor una raíz unitaria que se repite \(d\) veces: \(\left(1 - \sum_{i=1}^p \alpha_i L^i\right) = \left(1 - \sum_{i=1}^{p-d} \phi_i L^i\right) \cdot (1 - L^1)^d\). De esta forma ya podríamos definir nuestro modelo ARIMA(p’, q, d) como: \[ \left(\left(1 - \sum_{i=1}^{p'} \phi_i L^i\right) \cdot (1 - L^1)^d \right)X_t = \left(1 + \sum_{i=1}^q \theta_i L^i\right) \epsilon_t\\ con\ p' = p - d \] Si ahora hacemos \(p'=2, d=1, q=1\) nos queda ARIMA(2,1,1): \[ \left(\left(1 - \phi_1 L^1 - \phi_2 L^2\right) \cdot (1 - L^1)\right)X_t = \left(1 + \theta_1 L^1\right) \epsilon_t\ \\ X_t = \phi_1 X_{t-1} + \phi_2 X_{t-2} + X_{t-1} - \phi_1 X_{t-2} - \phi_2 X_{t-3} + \theta_1 \epsilon_{t-1} + \epsilon_t \\ X_t - X_{t-1} = \phi_1 \left(X_{t-1} - X_{t-2}\right)+ \phi_2 \left(X_{t-2} + X_{t-3}\right) + \theta_1\epsilon_{t-1} + \epsilon_t \\ \] Es importante señalar que el número de parámetros a ajustar para obtener el modelo es \(p + q\) (3 en nuestro caso). Para las series temporales con componente estacional se anaden los índices (P, Q, D) que cumplen la misma función que (p, q, d) para la parte no estacionaria pero sobre la componente estacional.

A continuación se presentan los residuos de un modelo ARMA estimado para la serie del D2ln(PIB). ¿Considera que dicho modelo es aceptable para realizar predicciones?. Justifique su respuesta.
Residuo

Residuo

Autocorrelación Residuo

Autocorrelación Residuo

Autocorrelación Parcial Residuo

Autocorrelación Parcial Residuo

Mirando los p-valores tanto para el test Box-Ljung como para el BDS podemos ver que todos los valores salvo el del test BDS para 3 dimensiones embebidas, son superiores a 0.05, luego no podemos rechazar la hipótesis nula en ninguno de los casos. De esta forma diremos que los residuos de nuestro modelo están independientemente distribuidos (test Ljung-Box) o indenpendiente e idénticamente distribuidos (test BDS) lo que concuerda con los gráficos de autocorrelación proporcionados donde para ningun orden (lag) encontramos valores elevados, y con la imagen de los propios residuos, donde vemos valores aleatorios distribuidos normalmente. Si los residuos se comportan como ruido blanco es muy probable que el modelo prediga con bastante solvencia el comportamiento de la serie, ya que debe estar recogiendo la parte explicable de su dinámica.

Ejercicio 10

De los diferentes indicadores adelantados para la predicción del PIB que se vieron en clase:

  • IndSintActEc : “I. Sintético Actividad (CVEC).Indice (2010=100)”

  • AfiliadoSS: “Afiliados Seguridad Social. Total. (Miles)”

  • ConsEnergEle: “Consumo energía eléctrica (CCal y Temp.). Millones de Kwh.”

  • IndSentEc: “I. Sent. Económico (1990 2017=100)”, Fuente Eurostat, promedio Indices Confianza Industrial, Servicios, Consumidores, Construción y Comercio

  • VentasGrdsEmpr: “Total Ventas grandes empresas.(CVEC).Indice (2013=100).”

  • ICC: Indice de Confianza del Consumidor

  • ComerMenor: Índice de Comercio Minorista

  • OCUPADOS: Número de OCUPADOS (EPA)

¿Cuál de todos ellos cree que es el más adecuado para adelantar el comportamiento de la tasa de crecimiento interanual del PIB?. ¿Se atrevería a adelantar cuál será la tasa de crecimiento interanual del PIB para el cuarto trimestre de este 2018?

En primer lugar debo indicar que no he encontrado en los ficheros trabajados en clase la variable ComerMenor y en ninguno de los sets de datos un índice de comercio minorista, así que he decidido sustituir este indicador por la serie temporal del IBEX35, ya que me parece que podía ser de interés.

El código empleado para contestar a las preguntas lo adjunto a continuación, con él se pueden obtener los valores y gráficos que utilizaré en adelante para justificar mis respuestas. Esencialmente lo que he hecho ha sido coger los indicadores mencionados y elaborar con ellos un modelo de regresión lineal de sus series, o de sus series diferenciadas, sobre la tasa de crecimienciento interanual del PIB. He estudiado además la cointegración entre las series, y probado incluso algunos modelos ARMAX, no obstante mi tiempo no es infinito, este procedimiento estaba recogido entre los procedimientos empleados en clase, y aunque puedan aparecer relaciones espúreas (todas ellas lo son ya que en el estudio de la cointegración nunca obtuve un estadístico superior a los valores tabulados para el test de cointegración de Phillips–Ouliaris, lo cual me resultó francamente desalentador), lo resultados parecen bastante buenos. El procedimiento puede reproducirse mediante la ejecución del siguiente código.

# Algunas librarÍas generales ----
library(dplyr)  # Filtrar y seleccionar variables
library(ggplot2) # Para realizar gr?ficos
library(tseries) # Adf.test, kpss.test, bds.test, get.hist.quote, portfolio.optim, surrogate, arma, garch
library(forecast) # interpolar, ndiffs para saber el n?mero de diferencias, adf, ggCcf(x,y): [xt+k, yt], ggAcf, autoarima
library(zoo) #diferentes procesos para series temporales
library(stargazer) # para tablas de resultados
library(tempdisagg) #ta agrega (reduce frecuencia), td desagrega (aumenta frecuencias)
library(pxR) # para leer datos desde archivos .px pc-axis del INE pxR::read.px() 
library(eurostat) #para leer datos de eurostat
library(data.table)
library(xts) #para series temporales grandes
library(dygraphs)
library(outliers) #outlier, rm.outlier, scores, chisq.out.test # para detectar outliers o datos an?malos ojo serie estacionaria
library(seasonal) #Para extraer cmponente estacional
library(tempdisagg) #ta agrega (reduce frecuencia), td desagrega (aumenta frecuencias)
library(urca)# Para test de ra?ces Unitarias
library(ecm) #para modelos de correci?n de error
library(lmtest) #varios test de modelos de regresi?n
library(sandwich) #estimaci?n robusta matriz de varianzas y covarianzas
library(rugarch) # estimaci?n de modelos GARCH


#_______________________________________________________________________________________
# 1 Carga los datos ----
# 1.0 Carga de Datos desde archivos RData ----

load(file="DATOCU.RData") # aquí he guardado las serie de Ocupados y 
                          # y de Tasa de Paro de la EPA

# 1.1. Carga de datos desde archivos (csv, excel, etc) ----

# Ejemplo: Series temporales de Coyuntura Econ?mica
# Bolet?n de Coyuntura Económica Semanal del Ministriode Economía 
# http://serviciosede.mineco.gob.es/indeco/
# tengo que descargar la hoja y desproteger la
# la hoja Cuadro2Mensual (pestaña Revisar. grupo Cambios, bot?n Desproteger hoja) y grabar

library(readxl)
IndicesMineco<-read_excel("indicador.xls",sheet = "Cuadro2 Mensual",skip = 2)
# Selecciono
datosMineco<-IndicesMineco[,c(1,2,3,4,7,8)]
datosMineco.desc<- names(datosMineco)
names(datosMineco)<-c("Periodos","IndSintActEc","AfiliadoSS","ConsEnergEle","IndSentEc", "VentasGrdsEmpr" )
str(datosMineco)
datosMineco$IndSintActEc<-as.numeric(datosMineco$IndSintActEc)
datosMineco$AfiliadoSS<-as.numeric(datosMineco$AfiliadoSS)
datosMineco$ConsEnergEle<-as.numeric(datosMineco$ConsEnergEle)
datosMineco$IndSentEc<-as.numeric(datosMineco$IndSentEc)
datosMineco$VentasGrdsEmpr<-as.numeric(datosMineco$VentasGrdsEmpr)

# Periodos : "Periodos"                                                  
# IndSintActEc : "I. Sintético Actividad (CVEC).Indice (2010=100)"           
# AfiliadoSS: "Afiliados Seguridad Social. Total. (Miles)"                
# ConsEnergEle: "Consumo energía el?ctrica (CCal y Temp.). Millones de Kwh."
# IndSentEc: "I. Sent. Económico (1990 2017=100)", Fuente Eurostat, promedio Indices Confianza Industrial, Servicios, Consumidores, Construci?n y Comercio
# VentasGrdsEmpr: "Total Ventas grandes empresas.(CVEC).Indice (2013=100)."

# Construcción del Indice Sintético de Actividad Económica
  # VGE: Ventas en grandes empresas 0,89
  # AFI: Afiliados a la Seguridad Social no agrarios 0,88
  # CEM: Consumo aparente de cemento 0,86
  # IMPNE: Importaciones no energ?ticas 0,67
  # ELE: Consumo de energ?a el?ctrica 0,62
  # IPI: Índice de producci?n industrial 0,60
  # ENT: Entrada de turistas 0,54
  # REP: Transporte RENFE pasajeros 0,54
  # AER: Tráfico aéreo de pasajeros 0,49
  # ICE: Indicador de Clima Econ?mico 0,42
  # REM: Transporte RENFE mercanc?as 0,21
  # FIN: Cr?dito a empresas y familias deflactado por IPC subyacente 0,17

# 1.2. Descarga de datos Online desde el INE ----
library(pxR) # para leer datos desde archivos .px pc-axis del INE pxR::read.px() 

#   References: http://www.ine.es/dyngs/INEbase/listaoperaciones.htm 
#   a) Series del PIB: Contabilidad Nacional TRIMESTRAL

#   INEbase < Economía < Cuentas Econ?micas < Contabilidad nacional trimestral (CNTR) de Espa?a <  Informaci?n detallada
#   Pestalla resultados (serie base 2010); Enlaces relacionados (serie bases anteriores)

url<-"http://www.ine.es/jaxiT3/files/t/es/px/28603.px?nocab=1"
datos<-as.data.frame(read.px(url), stringsAsFactor=FALSE)  #como un data.frame
table(datos$Agregados.macroeconómicos)
Pib<-datos%>%
  dplyr::filter(Niveles.y.tasas =="Dato base")%>%
  dplyr::filter(Agregados.macroeconómicos=="Producto interior bruto a precios de mercado")%>%
  dplyr::mutate(ANO=as.numeric(substr(Periodo,1,4)))%>%
  dplyr::mutate(Tipo.de.dato=ifelse(Tipo.de.dato=="Datos no ajustados de estacionalidad y calendario", "Bruto", "Desestacionalizado"))%>%
  dplyr::mutate(TRIM=as.numeric(substr(Periodo,6,7)))%>%
  dplyr::mutate(tiempo=ANO+TRIM*0.25-0.25)%>%
  dplyr::arrange(Tipo.de.dato,tiempo)%>%
  dplyr::select(Periodo, Tipo.de.dato, ANO, TRIM, tiempo, value)

# 1.3. Descarga de datos desde Eurostat ----
library(eurostat)
# Busco datos por tablas de Eurostat, ejemplo: PIB
query<-search_eurostat("GDP", type="table", fixed=T)
# View(query[1:5,])

# Busco los indicadores de Confianza
#http://appsso.eurostat.ec.europa.eu/nui/show.do?dataset=ei_bssi_m_r2&lang=en#

query<-search_eurostat("Sentiment", type="all")
# View(head(query,10))
temp.Eurostat<-get_eurostat("ei_bssi_m_r2",time_format = "num",type ="code", stringsAsFactors = FALSE,cache = FALSE) #,time_format = "raw",type ="code", stringsAsFactors = TRUE)
#Para ver que significan los campos
label_eurostat_vars(temp.Eurostat)
# aÑadir los labels
temp.Eurostat<-label_eurostat(temp.Eurostat, code=names(temp.Eurostat))
# Ahora tendr?a que seleccionar el dato de Espa?a, Francia, Italia, Reino Unido y Alemania, por ejemplo
# View(temp.Eurostat %>%
# distinct(indic, indic_code))

# BS-CCI-BAL: Construction confidence indicator
# BS-CSMCI-BAL: Consumer confidence indicator
# BS-ESI-I: Economic sentiment indicator
# BS-ICI-BAL: Industrial confidence indicator
# BS-RCI-BAL: Retail confidence indicator
# BS-SCI-BAL: Services Confidence Indicator

Sent.Eurostat<-temp.Eurostat %>%
  filter(geo_code %in% c("ES")) %>% 
  filter(indic_code=="BS-CSMCI-BAL") %>%
  filter(s_adj_code=="SA") %>%
  arrange(time_code)

# 1.4. Descarga de datos desde Yahoo Finanzas (librer?a quantmod)----
library(tseries) # get.hist.quote, portfolio.optim, surrogate, arma, garch
# descargamos la serie del IBEX35
# desde yahoo finanzas diarios, o desde Alpha Vantage (quantmod) datos minuto a minuto
Ibex <- get.hist.quote(instrument="^ibex", retclass="zoo")

# 2. Formateo de las series ----
# Extraigo de los datos cargados todas las series que voy a emplear y les asigno la misma frecuencia
# Creo dos versiones de la serie del PIB
datos<- Pib %>% 
  filter(Tipo.de.dato=="Bruto") %>% 
  arrange(tiempo) %>% 
  select(tiempo,value)

Pib.ts<-ts(data = datos$value, start = c(1995,1), frequency = 4)

datos<-Pib %>% 
  filter(Tipo.de.dato=="Desestacionalizado") %>% 
  arrange(tiempo) %>% 
  select(tiempo,value)

Pib.SA.ts<- ts(data = datos$value, start = c(1995,1), frequency = 4)
objetivo.SA<- datos$value[nrow(datos)-1]

# Creo la serie De Ocupados
Ocu.ts<-ts(DATOCU$OCU, start=1976.50, frequency = 4)

# Creo la serie De ICC
ICC.ts<-ts(Sent.Eurostat$values, start=min(Sent.Eurostat$time), frequency = 12)
ICC.ts<- ta(ICC.ts, conversion="last", to='quarterly', rm.na=TRUE)

# Creo la serie del IBEX
Ibex.zoo <- aggregate(Ibex$Close, as.yearqtr, mean)

# Creo las series Mineco
IndicesMineco.ts<- ts(datosMineco[,2:6], start=1995, frequency = 12)
IndicesMineco.ts<- ta(IndicesMineco.ts, conversion="average", to='quarterly')

# Juntamos las series (comando intersec, para eliminar valores posteriores al tercer trimestre de 2018)
All.ts<- ts.intersect(IndicesMineco.ts, Ocu.ts, ICC.ts, as.ts(Ibex.zoo), Pib.ts)
colnames(All.ts)<- c("IndSintActEc", "AfiliadoSS", "ConsEnergEle", "IndSentEc", "VentasGrdsEmpr",
                     "Ocu.ts", "ICC.ts", "Ibex.ts", "Pib.ts")

# 3. Gráfico de las series ----
# plot(serie_temporal)  en realidad llama al procedimiento plot.ts, plot.zoo, plot.xts
colfunc<-colorRampPalette(c("firebrick","goldenrod","darkolivegreen3","cornflowerblue", "darkviolet"))
plot(scale(All.ts), plot.type = "single", col=colfunc(ncol(All.ts)))
legend('topleft','groups',colnames(All.ts), lty=1, cex=0.65, lwd=2.5, col=colfunc(ncol(All.ts)))

library(ggplot2)
aux<- as.data.frame(scale(All.ts))[index(All.ts)>2000,]
aux$id = 1:nrow(aux)
aux<- melt(aux, id='id')
ggplot() + geom_line(data = aux, aes(x = id, y = value, color = variable, group = variable), size = 1) +
  labs(title="Todas las series",
       subtitle = "2000Q1 - 2018Q3",
       y= "Valores escalados",
       x= "trimestre")

# 4. Análisis de datos Missing e imputación ----
# Cojo todas las series desde el primer dato no missing, salvo la del IBEX que la cojo un no missing posterior, ya que aún hay
# regiones de timepo amplias sin datos tras el primer missing. Los datos missing intermedios los interpolo.
# Además busco outliers, valores a más de 3 desviaciones típicas de la media

# INDICE DE ACTIVIDAD ECONÓMICA
IndSintActEc.ts<- IndicesMineco.ts[,1]
table(is.na(IndSintActEc.ts)) # No hay missing
which(abs(scores(IndSintActEc.ts, type="z"))>3) # No hay outliers

# AFILIADOS A LA SS
AfiliadoSS.ts<- IndicesMineco.ts[,2]
table(is.na(AfiliadoSS.ts)) # Eliminamos los missing
AfiliadoSS.ts<- window(AfiliadoSS.ts, start=index(AfiliadoSS.ts)[which.min(is.na(AfiliadoSS.ts))])
which(abs(scores(AfiliadoSS.ts, type="z"))>3) # No hay outliers

# CONSUMO ENERGÉTICO
ConsEnergEle.ts<- IndicesMineco.ts[,3]
table(is.na(ConsEnergEle.ts)) # Eliminamos los missing
ConsEnergEle.ts<- window(ConsEnergEle.ts, start=index(ConsEnergEle.ts)[which.min(is.na(ConsEnergEle.ts))])
which(abs(scores(ConsEnergEle.ts, type="z"))>3)

# INDICADOR DE SENTIMIENTO ECONÓMICO
IndSentEc.ts<- IndicesMineco.ts[,4]
table(is.na(IndSentEc.ts)) # No hay missing
which(abs(scores(IndSentEc.t, type="z"))>3) # No hay outliers

# VENTAS DE GRANDES EMPRESAS
VentasGrdsEmpr.ts<- IndicesMineco.ts[,5]
table(is.na(VentasGrdsEmpr.ts)) # Eliminamos los missing
VentasGrdsEmpr.ts<- window(VentasGrdsEmpr.ts, start=index(VentasGrdsEmpr.ts)[which.min(is.na(VentasGrdsEmpr.ts))])
which(abs(scores(VentasGrdsEmpr.ts, type="z"))>3) # No hay outliers

# PERSONAS OCUPADAS
table(is.na(Ocu.ts)) # No hay missing
which(abs(scores(Ocu.ts, type="z"))>3) # No hay outliers

# INDICE DE CONFIANZA DEL CONSUMIDOR
table(is.na(ICC.ts)) # No hay missing
which(abs(scores(ICC.ts, type="z"))>3) # No hay outliers

# IBEX35
Ibex.ts<- as.ts(Ibex.zoo)
table(is.na(Ibex.ts)) # Eliminamos los missing del principio
Ibex.ts<- window(Ibex.ts, start=c(2005,2)) # Inputamos missing
plot(Ibex.ts, lwd=2)
y1<- zoo::na.approx(Ibex.ts)
lines(y1, col="blue",type="p", pch=15, cex=1)
y2<-zoo::na.locf(Ibex.zoo)
lines(y2, col="darkgreen",type="p", pch=17, cex=0.9)
y<- Ibex.ts
y[1]<- zoo::na.approx(Ibex.ts)[1]
y3<-zoo::na.StructTS(y)
lines(y3, col="red",type="p", pch=19, cex=0.8) # Select this
Ibex.ts<- y3
which(abs(scores(Ibex.ts, type="z"))>3) # No hay outliers

# PIB
table(is.na(Pib.ts)) #Elliminamos los missing
which(abs(scores(Pib.ts, type="z"))>3) # No hay outliers


# 5. Desestacionalización ----
# Tendencia + Ciclo + Estacionalidad + SemanaSanta y Efecto Calendario + Errático
library(seasonal)
y<-IndSintActEc.ts
y_trend1<-decompose(y) #se limita a hacer una media movil
y_trend2<-seas(y) # método X13-Seats
plot(y)
lines(y_trend1$trend, col=2)
lines(trend(y_trend2), col=4)# Esta no presenta componente estacional

y<-AfiliadoSS.ts
y_trend1<-decompose(y) #se limita a hacer una media movil
y_trend2<-seas(y) # método X13-Seats
plot(y)
lines(y_trend1$trend, col=2)
lines(trend(y_trend2), col=4)# Esta presenta componente estacional
AfiliadoSS.ts<- trend(y_trend2)


y<-ConsEnergEle.ts
y_trend1<-decompose(y) #se limita a hacer una media movil
y_trend2<-seas(y) # método X13-Seats
plot(y)
lines(y_trend1$trend, col=2)
lines(trend(y_trend2), col=4)# Esta presenta componente estacional
ConsEnergEle.ts<- y_trend1$trend
ConsEnergEle.ts[is.na(ConsEnergEle.ts)]<- trend(y_trend2)[is.na(ConsEnergEle.ts)]

y<-IndSentEc.ts
y_trend1<-decompose(y) #se limita a hacer una media movil
y_trend2<-seas(y) # método X13-Seats
plot(y)
lines(y_trend1$trend, col=2)
lines(trend(y_trend2), col=4)# Esta no presenta componente estacional

y<-VentasGrdsEmpr.ts
y_trend1<-decompose(y) #se limita a hacer una media movil
y_trend2<-seas(y) # método X13-Seats
plot(y)
lines(y_trend1$trend, col=2)
lines(trend(y_trend2), col=4)# Esta no presenta componente estacional

y<-Ocu.ts
y_trend1<-decompose(y) #se limita a hacer una media movil
y_trend2<-seas(y) # método X13-Seats
plot(y)
lines(y_trend1$trend, col=2)
lines(trend(y_trend2), col=4)# Esta presenta componente estacional
Ocu.ts<- trend(y_trend2)

y<-ICC.ts
y_trend1<-decompose(y) #se limita a hacer una media movil
y_trend2<-seas(y) # método X13-Seats
plot(y)
lines(y_trend1$trend, col=2)
lines(trend(y_trend2), col=4)# Esta presenta componente estacional
ICC.ts<- trend(y_trend2)

y<-Ibex.ts
y_trend1<-decompose(y) #se limita a hacer una media movil
y_trend2<-seas(y) # método X13-Seats
plot(y)
lines(y_trend1$trend, col=2)
lines(trend(y_trend2), col=4)# Esta no presenta componente estacional

y<-Pib.ts
y_trend1<-decompose(y) #se limita a hacer una media movil
y_trend2<-seas(y) # método X13-Seats
plot(y)
lines(y_trend1$trend, col=2)
lines(trend(y_trend2), col=4)
lines(Pib.SA.ts, col=3)# Esta si presenta componente estacional
Pib.ts<-Pib.SA.ts

# 6. Test de Raices Unitarias y cointegración ----
# Voy a crear dos series independientes
# Tienen algo que ver la una con la otra?
cointegradas<- function(seriex, seriey){
  # 1 Estacionariedad
  series<-ts.intersect(seriey, seriex)
  y = series[,1]
  x = series[,2]
  
  ggtsdisplay(x)
  ggtsdisplay(y)
  
  # test de Raices Unitarias
  library(tseries)
  print(adf.test(x))
  print(adf.test(y))
  
  # Opción uno Diferencio las series y compruebo que el modelo siga siendo válido
  #construyo las series diferenciadas dy dx
  cat("Serie x diferenciada", ndiffs(x), "veces\n")
  cat("Serie y diferenciada", ndiffs(y), "veces\n")
  dy <- diff(y, ndiffs(y))
  dx <- diff(x, ndiffs(x))
  series.dif<- ts.intersect(dy,dx)
  dy<- series.dif[,1]
  dx<- series.dif[,2]
  
  # Construyo las regresiones entre las series, y entre las series diferenciadas
  coint.reg <- lm(y ~ x)
  summary(coint.reg)
  dif.reg<-lm(dy~dx)
  print(summary(dif.reg))
  library(stargazer)
  stargazer(coint.reg, dif.reg, type="text")
  
  # Lo ploteamos
  aux<- data.frame(real=y, ajustado=coint.reg$fitted)
  aux$id = 1:nrow(aux)
  aux<- melt(aux, id='id')
  print(ggplot() + geom_line(data = aux, aes(x = id, y = value, color = variable, group = variable), size = 1))
  
  #Test cointegración Engle-Granger
  # Primera etapa
  # extraigo el error
  coint.err<-residuals(coint.reg)
  plot(coint.err, type="l")
  
  # A menudo la regresión sobre las series diferenciadas SI pasan los test de cointegración !!!
  # segunda etapa: compruebo si el error es estacionario
  cat("Augmented Dickey Fuller test p-value:", punitroot(adf.test(coint.err)$statistic, trend="nc"), "\n")
  
  # Phillips–Ouliaris Cointegration Test HO: NO Cointregration 
  # null hypothesis that x is not cointegrated.
  library(urca)
  # null hypothesis that x is not cointegrated.
  tru<-ur.df(coint.err, type = c("drift"), lags = 12, 
             selectlags = c("BIC"))
  print(summary(tru))
  cat("Augmented Dickey Fuller test p-value:", punitroot(tru@teststat[1], trend="c"), "\n")
  capocointest<-ca.po(series, demean=c("none"), lag=c("long"))
  summary(capocointest)
}

seriey<- Pib.ts
seriex<- IndSintActEc.ts
cointegradas(seriex, seriey)
seriex<- AfiliadoSS.ts
cointegradas(seriex, seriey)
seriex<- IndSentEc.ts
cointegradas(seriex, seriey)
seriex<- ConsEnergEle.ts
cointegradas(seriex, seriey)
seriex<- VentasGrdsEmpr.ts
cointegradas(seriex, seriey)
seriex<- Ocu.ts
cointegradas(seriex, seriey)
seriex<- ICC.ts
cointegradas(seriex, seriey)
seriex<- Ibex.ts
cointegradas(seriex, seriey)


# 7. Modelos de seres temporales Estacionarias ----
# Correlación contemporanea e indicadores adelantados
# Tasa de crecimiento interanual del Pib
obj.comp<- Pib.ts/stats::lag(Pib.ts, -4) -1
objetivo<- as.data.frame(window(obj.comp, start=c(2018,4)))[1,]
seriey<- window(obj.comp, end=c(2018,3))
lin.model<- function(seriey, seriex, nombre){
  ggCcf(seriex, seriey)
  crosscorr<- Ccf(seriex, seriey)
  Ccf.df<- data.frame(values=crosscorr$acf, row.names=crosscorr$lag)
  n_lags<- as.numeric(row.names(Ccf.df)[which.max(Ccf.df$values)])
  
  # Test de Causalidad de Granger 
  library(lmtest)
  grangertest(seriex, seriey, order=2)
  grangertest(seriey, seriex, order=2)
  
  if(n_lags>0){
    cat(nombre, " es un indicador atrasado ", n_lags, " periodos\n")
  }else if(n_lags<0){
    cat(nombre, " es un indicador adelantado ", n_lags, " periodos\n")
  }else{
    cat(nombre, " es un indicador contemporáneo ", n_lags, " periodos\n")
  }
  seriex2<-stats::lag(seriex, n_lags)
  series<-ts.intersect(seriey,seriex2)
  modelo.lm<-lm(seriey ~ seriex2, data=series)  #ojo, na.action=NULL para las series
  cat("El modelo asociado a la serie lagueada presenta un R2 ajustado de: ", summary(modelo.lm)$adj.r.squared, "\n")
  to.pred<- as.data.frame(window(seriex2, start=c(2018,4)))[1,]
  predicho<- as.numeric(modelo.lm$coefficients[1]) + as.numeric(modelo.lm$coefficients[2])*to.pred
}

# Adaptamos la seriex para que se parezca a la seriey
seriex<- IndSintActEc.ts
plot(scale(seriey), type="l")
lines(scale(seriex), type="l", col=2)
lines(scale(diff(seriex)), type="l", col=4)
# Usamos la seriex que nos convenga
seriex<- diff(seriex)
predicho<- lin.model(seriey, seriex, "IndSintActEc.ts")
cat("Error: ", abs(predicho - objetivo)*100/objetivo, "%\n")

seriex<- AfiliadoSS.ts
plot(scale(seriey), type="l")
lines(scale(seriex), type="l", col=2)
lines(scale(diff(seriex)), type="l", col=4)
# Usamos la seriex que nos convenga
seriex<- diff(seriex)
predicho<- lin.model(seriey, seriex, "AfiliadoSS.ts")
cat("Error: ", abs(predicho - objetivo)*100/objetivo, "%\n")

seriex<- IndSentEc.ts
plot(scale(seriey), type="l")
lines(scale(seriex), type="l", col=2)
lines(scale(diff(seriex)), type="l", col=4)
# Usamos la seriex que nos convenga
predicho<- lin.model(seriey, seriex, "IndSentEc.ts")
cat("Error: ", abs(predicho - objetivo)*100/objetivo, "%\n")

seriex<- ConsEnergEle.ts
plot(scale(seriey), type="l")
lines(scale(seriex), type="l", col=2)
lines(scale(diff(seriex)), type="l", col=4)
# Usamos la seriex que nos convenga
seriex<- diff(seriex)
predicho<- lin.model(seriey, seriex, "ConsEnergEle.ts")
cat("Error: ", abs(predicho - objetivo)*100/objetivo, "%\n")

seriex<- VentasGrdsEmpr.ts
plot(scale(seriey), type="l")
lines(scale(seriex), type="l", col=2)
lines(scale(diff(seriex)), type="l", col=4)
# Usamos la seriex que nos convenga
seriex<- diff(seriex)
predicho<- lin.model(seriey, seriex, "VentasGrdsEmpr.ts")
cat("Error: ", abs(predicho - objetivo)*100/objetivo, "%\n")

seriex<- Ocu.ts
plot(scale(seriey), type="l")
lines(scale(seriex), type="l", col=2)
lines(scale(diff(seriex)), type="l", col=4)
# Usamos la seriex que nos convenga
seriex<- diff(seriex)
predicho<- lin.model(seriey, seriex, "Ocu.ts")
cat("Error: ", abs(predicho - objetivo)*100/objetivo, "%\n")

seriex<- ICC.ts
plot(scale(seriey), type="l")
lines(scale(seriex), type="l", col=2)
lines(scale(diff(seriex)), type="l", col=4)
# Usamos la seriex que nos convenga
seriex<- diff(seriex)
predicho<- lin.model(seriey, seriex, "ICC.ts")
cat("Error: ", abs(predicho - objetivo)*100/objetivo, "%\n")

seriex<- Ibex.ts
plot(scale(seriey), type="l")
lines(scale(seriex), type="l", col=2)
lines(scale(diff(seriex)), type="l", col=4)
# Usamos la seriex que nos convenga
seriex<- diff(seriex)
predicho<- lin.model(seriey, seriex, "Ibex.ts")
cat("Error: ", abs(predicho - objetivo)*100/objetivo, "%\n")
# 8. Los mejores ----
lin.model2<- function(seriey, seriex, nombre){
  ggCcf(seriex, seriey)
  crosscorr<- Ccf(seriex, seriey)
  Ccf.df<- data.frame(values=crosscorr$acf, row.names=crosscorr$lag)
  n_lags<- as.numeric(row.names(Ccf.df)[which.max(Ccf.df$values)])
  
  # Test de Causalidad de Granger 
  library(lmtest)
  grangertest(seriex, seriey, order=2)
  grangertest(seriey, seriex, order=2)
  
  if(n_lags>0){
    cat(nombre, " es un indicador atrasado ", n_lags, " periodos\n")
  }else if(n_lags<0){
    cat(nombre, " es un indicador adelantado ", n_lags, " periodos\n")
  }else{
    cat(nombre, " es un indicador contemporáneo ", n_lags, " periodos\n")
  }
  seriex2<-stats::lag(seriex, n_lags)
  series<-ts.intersect(seriey,seriex2)
  modelo.lm<-lm(seriey ~ seriex2, data=series)  #ojo, na.action=NULL para las series
  cat("El modelo asociado a la serie lagueada presenta un R2 ajustado de: ", summary(modelo.lm)$adj.r.squared, "\n")
  to.pred<- data.frame(row.names=index(seriex2), seriex2)["2019.25",]
  as.numeric(modelo.lm$coefficients[1]) + as.numeric(modelo.lm$coefficients[2])*to.pred
}
seriex<- AfiliadoSS.ts
seriex<- diff(seriex)
pred1<- lin.model2(obj.comp, seriex, "AfiliadosSS")
seriex<- ICC.ts
seriex<- diff(seriex)
pred2<- lin.model2(obj.comp, seriex, "AfiliadosSS")
aux<- stats::lag(Pib.ts, -4)
pred<- (1 + (pred1+pred2)/2)*data.frame(row.names=index(aux), aux)["2019.25",]
cat("El valor predicho para el segundo trimestre de 2019 del Pib es: ", pred)

A la vista de los resultados coincideremos en que todos los indicadores son indicadores que anticipan el comporamiento del PIB, incluido el indicador de Sentimiento Económico, lo cual contradice el comportamiento observado en la gráficas. Mirando ya los resultados, dado que actualmente los datos del PIB contienen su valor en el 4º trimestre de 2018, he podido ver que error se comete al predecir este valor con los distintos modelos, considerando solo esta cifra podría uno puede pensar que el Índice de Confianza del Consumidor nos da el mejor modelo, ya que sólo se desvió un 7% del valor real, sin embargo, si miramos también el \(R^2\), debemos quedarnos con el modelo asociado a las diferencias de la serie de Número de Afiliados a la Seguridad Social, lagueadas 2 periodos (6 meses). Este modelo presenta un \(R^2\) de 0.896 lo que resulta un valor notable, y su predicción para el 4º trimestre de 2018 tan solo se desvió un 26.456% del valor real, de esta forma el valor que habría sugerido para Tasa de Crecimiento Interanual del PIB sería de un: 0.0294% mientras que la real fue 0.0233%.

Por no escaquearme del juego de las predicciones diré que para el 2º trimestre de 2019 espero una tasa del: 0.0323%, luego el Pib valdrá: 111.384, lo cual implica que aumentará más de 1.5 puntos respecto a este trimestre, lo que supone la mayor subida de la que se tienen datos, ergo, es un dato seguramente muy optimista. Si promedio el valor que nos da el modelo con mejor \(R^2\) (el asociado a los Afiliados a la Seguridad Social) con el que mejor predijo el valor para el último trimestre de 2018 (el asociado al Índice de Confianza de los Consumidores) tenemos una tasa de crecimiento del 0.0269% y un valor del PIB de 110.81 puntos, que ya suena más realista.

Riesgo de Crédito y Scoring

Ejercicio 11

Defina qué es el concepto de riesgo de crédito y qué es un tarjeta de puntuación o Scoring

El riesgo de crédito se engloba dentro del riesgo financiero que está a su vez englobado dentro de los riesgos económicos. El riesgo de crédito contempla el riesgo asociado al incumplimiento de sus obligaciones por alguna de las partes involucradas en una operación financiera, así como del riesgo de que la propia operación sea fraudulenta. Una tarjeta de Scoring califica a los posibles deudores, es decir, les asigna una puntuación en función del riesgo existente a que cometan impago.

Ejercicio 12

Suponga que está realizando un modelo de Scoring de riesgo de impago de créditos, y que dispone de la siguiente información sobre el porcentaje de créditos impagados (m2av) según el tiempo de vigencia del préstamo desde su concesión inicial ¿qué horizonte temporal utilizaría para definir la variable objetivo?. Razone su respuesta
Créditos impagados



A la vista del gráfico sabemos que la mayor parte de los créditos impagados se producen entre 24 y 36 meses tras la concesión del préstamo. Por lo tanto el periodo de estudio que emplearía sería aquel que contuviera mayor cantidad de información de préstamos en este periodo de vigencia. De esta forma aumento mis probabilidades de conseguir información sobre las características de un impago. Una forma aún más precisa de hacerlo sería asignar a cada periodo un peso proporcional a la curva observada y así asignar una puntuación a cada préstamo en función del periodo en que se encuentre, y luego hacer una ventana temporal que sume la mayor puntuación posible, que contenga la mayor cantidad de observaciones de préstamos entre 20 y 40 meses después de su concesión.

Ejercicio 13

Una de las fases más importantes en el proceso de la creación de la tarjeta de puntuación es el de la primera selección de las potenciales variables explicativas que podrían incluirse en el modelo de regresión. Por favor explique al menos dos estadísticos que podrían utilizarse para esta tarea de selección inicial, que significan y cómo deben usarse

Para seleccionar variables podemos emplear diferentes estadísticos, que determinan la capacidad de predicción de la variable sobre una variable objetivo:

  • Valor de Información \((IV)\): Se emplea con variables cualitativas (o cuantitativas discretizadas) y nos permite ver la información aportada por cada una de las categorías para la discriminación de la variable objetivo, podemos además categorizar las variables según si IV:

    • IV \(\lt\) 0.02 –> Variables no predictivas
    • 0.02 \(\leq\) IV \(\lt\) 0.1 –> Variables de influencia débil
    • 0.1 \(\leq\) IV \(\lt\) 0.3 –> Variables de influencia media
    • 0.3 \(\leq\) IV –> Variables de influencia fuerte

    Matemáticamente se define según:

\[ IV = \sum_{i=1}^{nº\ valores}WE_i \cdot \left( \frac{G_i \cdot B - B_i \cdot G}{B \cdot G} \right)\\ \text{siendo: }\ WE_i = \ln \left( \frac{G_i \cdot B}{B_i \cdot G} \right)\ \text{, } G = goods \text{ y }B = bads \]

  • Coeficiente de Gini \((Gini)\): Se desarrollo como un índice para determinar la desigualdad económica, sin embargo, a día de hoy su uso también se ha ampliado. Igual que el estadístico anterior se emplea únicamente con variables cualitativas (o cuantitativas discretizadas) y también nos permite categorizar las variables según:
    • G \(\lt\) 5%: No se recomienda el uso de la variable
    • 5% \(\leq\) G \(\lt\) 15%: Se puede emplear la variable
    • 15% \(\leq\) G: Se recomienda el uso de la variable
    Matemáticamente se define según:

\[ Gini = 1 - \frac{2 \cdot \sum_{i=2}^m \left( G_i \cdot \sum_{j=1}^{i-1} B_j \right) + \sum_{k=1}^{m} \left(G_k B_k\right)}{G \cdot B}\\ \text{siendo: } G = goods \text{ y }B = bads \]

  • V de Cramer \((\Phi_C)\): Al igual que los anteriores se emplea con variables cualitativas (o cuantitativas discretizadas) y también nos permite categorizar las variables según:
    • \(\Phi_C\) \(\lt\) 0.2: No hay relación entre las variables
    • 0.2 \(\leq\) \(\Phi_C\) \(\lt\) 0.6: Existe una relación moderada entre las variables
    • 0.6 \(\leq\) \(\Phi_C\): Existe una relación fuerte entre las variables
    Matemáticamente se define, si llamamos \(k\) al número de categorías de la variable objetivo, y \(l\) al de la variable estudiada, como:

\[ \Phi_C = \sqrt{\frac{\chi^2}{N \cdot (\min(k, l) -1)}}\\ \text{siendo: } \chi^2 = \sum_{i=1}^k \sum_{j=1}^l \frac{\left( n_{ij} - \tfrac{n_i n_j}{n} \right)^2}{\frac{n_i n_j}{n}} \]

Ejercicio 14

Una vez depurada la base de datos, analizados atípicos e imputados, en su caso, los valores perdidos, una vez definida la variable objetivo, dividida la base de datos en muestra de entrenamiento y de validación, y una vez elegidas las variables explicativas que se incluirán en el modelo, procede la estimación de un modelo de probabilidad. Defina y explique las características de un modelo de probabilidad o de variable dependiente binaria, y las limitaciones de este tipo de modelos para realizar la diagnosis de sus residuos y la interpretación de sus coeficientes

Un modelo de probabilidad o variable dependiente binaria es un modelo que diseñado para predecir una variable binaria a partir otras variables. El modelo es por lo tanto una función de enlace (en general una distribución de probabilidad) \(G(\overrightarrow{X})\) sobre el espacio de las variables, el valor de esta función puede considerarse una variable latente, y la variable objetivo binaria simplemente una función sobre los valores de la función de enlace que toma los valores uno o cero, tenemos varios modelos según su función de enlace: - Modelo Lineal de Probabilidad: Distribución Uniforme - Modelo Probit: Distribución normal - Modelo Logit: Distribución logística - Modelo Gompit: Valor extremo

En un modelo logit o de regresión logística, que es el más común, la distribución de probabilidad es: \[ P(y=1 | \vec{X}) = \frac{\exp \left(\beta_0 + \vec{\beta}\cdot \vec{X} \right)}{1 + \exp \left(\beta_0 + \vec{\beta}\cdot \vec{X} \right)}\\ P(y=0 | \vec{X}) = 1 - P(y=1 | \vec{X}) = \frac{1}{1 + \exp \left(\beta_0 + \vec{\beta}\cdot \vec{X} \right)}\\ odds = \frac{P(y=1 | \vec{X})}{P(y=0 | \vec{X})} = \exp \left(\beta_0 + \vec{\beta}\cdot \vec{X} \right)\\ Logit(P) = \ln(odds) = \beta_0 + \vec{\beta}\cdot \vec{X} \longrightarrow Lineal \] En este caso los coeficientes admiten una interpretación directa (de ahí que se haya extendido su uso), ya que su exponencial representa el odds ratio asociado al aumento de la variable una unidad. En modelos con distribuciones diferentes la interpretación de los parámetros implica mayor complejidad. A la hora de estimar los parámetros no podemos emplear mínimos cuadrados, o alguna técnica analítica similar, sino que tenemos que recurrir al argumento de máxima verosimilitud que implica asumir que los datos que hemos tomado son los que son porque este evento era el más probable. De esta forma los parámetros se aproximan mediante modelos iterativos y esto implica que solo conozcamos el residuo de los modelos de forma aproximada.

Ejercicio 15

Relacionado con el punto anterior, comente qué es el porcentaje de pronósticos correctos, y cómo puede utilizarse el umbral para pronosticar si un préstamo será impagado o no.

Una vez construido el modelo hay varios valores que nos serán de interés:

  • Tasa de aciertos: \(\frac{VP + VN}{TODOS}\)

  • Tasa de fallos: \(\frac{FP + FN}{TODOS}\)

  • Sensibilidad (positivos acertados entre totales): \(\frac{VP}{VP + FN}\)

  • Especificidad (negativos acertados entre totales): \(\frac{VN}{VN + FP}\)

  • Valor predictivo positivo (positivos acertados entre pronosticados): \(\frac{VP}{VP + FP}\)

  • Valor predictivo negativo (negativos acertados entre pronosticados): \(\frac{VN}{VN + FN}\)

Dado que lo que tenemos es esa variable latente que asigna a cada observación, en función de los valores de algunas variables, un valor entre 0 y 1, podemos elegir cualquier punto de corte a partir del cual valores mayores sean uno y menores cero. A esto lo llamamos umbral, los prétamos que lo superen se calificarán como impagados. No obstante este umbral lo definimos nosotros, con lo que dependiendo el valor que elijamos tendremos mayor o menor número de préstamos pronosticados como impagados. Para elegir el punto de corte podemos también emplear varios criterios diferentes, en función de qué prioricemos:

  • El punto que maximiza la tasa de aciertos

  • El punto que maximiza el índice de Youden: \(Youden = sensibilidad + especificidad -1\)

  • El punto que, fijada una especifidad mínima, maximiza la sensibilidad

  • El punto que, fijada una sensibilidad mínima, maximiza la especificidad

Explique brevemente qué es la curva ROC y cómo debe interpretarse, ¿sobre que muestra cree que es mejor realizar la diagnosis utilizando esta curva, con la muestra de entrenamiento o con la de validación?. Justifique su respuesta.

Para conocer la bondad de nuestro modelo se emplea la curva ROC (acónimo de Receiver Operating Characteristic) que es una representación de la sensibilidad del modelo frente a uno menos la especificidad. Esto implica que cuanto más próxima esté a la esquina superior izquierda del gráfico (sensibilidad = 1, especificidad = 1) mejor será el modelo. Normalmente se emplea el área bajo la curva como un estimador de la bondad del modelo, también podríamos elegir como punto de corte el punto de la curva que se encontrase a menor distancia del vértice antes señalado.

Ejercicio 16

Explique brevemente qué es “la inferencia de denegados”, cuándo conviene realizarla y por qué

La inferencia de denegados es una técnica que se solventar el hecho de que la muestra empleada recoge únicamente los prétamos otorgados, y no aquellos que se rechazaron, es decir, de los préstamos que no otorgamos no sabemos si serían impagados o no, se produce así un sesgo derivado de la selección de la muestra que puede conducir al modelo a infravalorar la probabilidad de impago. Para salvar este sesgo se emplea la inferencia de denegados que se aplica fundamentalmente mediante cuatro técnicas, asignar todos los préstamos no otorgados a la categoría de impagados, asignar a los denegados la categoría que defina el modelo creado, estableciendo un ‘’cut-off’’, añadir una variable con la probabilidad asignada con el modelo generado con los préstamos aceptados multiplicada por un ratio definido por nosotros y desdoblar cada denegado en 2, uno con la probabilidad de impago y como si hubiera impagado, otro con la probabilidad de pago y como si hubiera pagado, por último esta la opción de asignarles un valor aleatorio pocedente de la distribución que presentan los prestamos aceptados que obtienen la misma probabilidad con el modelo de aceptados.

Ejercicio 17

Construya una tarjeta de puntuación de riesgos de impago con el conjunto de datos del fichero “datriesgos.RData” siguiendo la metodología propuesta en clase. El fichero contiene:
  • Variable objetivo:
    • Impago (1 SI impago, 0 NO impago)
  • Variables explicativas:
    • loan_amnt (Importe del pretamo)
    • int_rate (Interes aplicado)
    • calidad_crediticia (A - G siendo A menor riesgo)
    • emp_length (antiguedad del empleo)
    • home_ownership (1 propietario, 0 NO propietario)
    • annual_inc (ingresos anuales del cliente)
    • age (edad)

En este ejercicio he puesto bastante esfuerzo así que voy a ir comentando el código poco a poco. En primer lugar, como siempre, he cargado las librerías y los datos, le he echado un ojo, y he definido algunas funciones que usaré de forma recurrente. Conviene señalar que es más que recomendable leerse el código ya que está llenos de apuntes sobre los resultados.

# Descripción: Cálculo de una scorecard usando el paquete de R "scorecard"
setwd("~/Desktop/Santi/master_BD/temporal_analisis/entrega/")
# 0. Cargo librerias ----
# Traditional Credit Scoring Using Logistic Regression
library(stargazer)
library(mfx) # para efectos marginales
library(scorecard)
library(outliers)
library(gmodels) #CrossTable()
library(car)
library(ggplot2)
library(gridExtra)
library(onewaytests) #para el test welch
library(caret)
library(questionr)

## 0.1 Funciones útiles ----
Vcramer<-function(v,target){
  if (is.numeric(v)){
    v<-cut(v,5)
  }
  if (is.numeric(target)){
    target<-cut(target,5)
  }
  cramer.v(table(v,target))
}

graficoVcramer<-function(matriz, target, title){
  salidaVcramer<-sapply(matriz,function(x) Vcramer(x,target))
  barplot(sort(salidaVcramer,decreasing =T),las=2,main=title, horiz=TRUE,
          cex.axis=0.65, xlim=c(0,1))
}

myDistPlot<- function(datos, var1, var2, title){
  ggplot(datos, aes_string(x=var1, fill=var2)) +
    geom_density(alpha=.3, show.legend = FALSE) + ggtitle(title)
}

myBoxPlot<- function(datos, var1, var2, varbin, title, leg){
  ggplot(datos, aes_string(x=var1, y=var2, fill=varbin)) + geom_boxplot(alpha=.3, show.legend = leg) + 
    stat_summary(fun.y=mean, colour="darkolivegreen", geom="point", shape=18, size=3,
                 position=position_dodge(width=0.75), show_guide = FALSE) + 
    ggtitle(title)
}

# Rellena los missings en var1 con valores de var1 parecidos a entradas con valores similares de var2
rellenaMissings<- function(dat, var1, var2, no.lower){
  dat$var2.fac<- cut(dat[[var2]], 5)
  missings<- dat[is.na(dat[[var1]]),c(var1, var2, "var2.fac")]
  not.missings<- dat[!is.na(dat[[var1]]), c(var1, var2, "var2.fac")]
  inputator<- as.data.frame(aggregate(not.missings[,1], list(not.missings$var2.fac), mean))
  inputator$sd<- aggregate(not.missings[,1], list(not.missings$var2.fac), sd)$x
  set.seed(12345)
  if(no.lower){
    aleat<- runif(nrow(missings), 0, 1)
  }else{
    aleat<- runif(nrow(missings), -1, 1)
  }
  for(i in seq(1,nrow(missings))){
    miss.group<- missings$var2.fac[i]
    al<- aleat[i]
    val<- inputator[inputator$Group.1 == miss.group, 2] + al * inputator[inputator$Group.1 == miss.group, 3]
    dat[rownames(missings)[i],var1]<- val
  }
  return(dat[, -ncol(dat)])
}

# 1.Cargo los datos ----
load("datriesgos.RData")
datos<- datriesgos
#Visualización
datos

Como se puede observar se trata de un set con 8 variables y 29092 observaciones, luego la proporción nos va a permitir tener un modelo con un número de parámetros razonable (por el número de variables) y bien ajustado (por el número de observaciones).

A continuación voy a buscar outliers o datos que me parezcan erróneos, reemplazarlos por NAs y luego inputar valores en los datos missing.

## 1.1 Hago factores ----
datos$impago<- as.factor(datos$impago)
# summary
summary(datos)
##  impago      loan_amnt        int_rate     Calidad_Crediticia
##  0:25865   Min.   :  500   Min.   : 5.42   A:9649            
##  1: 3227   1st Qu.: 5000   1st Qu.: 7.90   B:9329            
##            Median : 8000   Median :10.99   C:5748            
##            Mean   : 9594   Mean   :11.00   D:3231            
##            3rd Qu.:12250   3rd Qu.:13.47   E: 868            
##            Max.   :35000   Max.   :23.22   F: 211            
##                            NA's   :2776    G:  56            
##    emp_length      home_ownership    annual_inc           age       
##  Min.   : 0.000   MORTGAGE:12002   Min.   :   4000   Min.   : 20.0  
##  1st Qu.: 2.000   OTHER   :   97   1st Qu.:  40000   1st Qu.: 23.0  
##  Median : 4.000   OWN     : 2301   Median :  56424   Median : 26.0  
##  Mean   : 6.145   RENT    :14692   Mean   :  67169   Mean   : 27.7  
##  3rd Qu.: 8.000                    3rd Qu.:  80000   3rd Qu.: 30.0  
##  Max.   :62.000                    Max.   :6000000   Max.   :144.0  
##  NA's   :809

## 1.2 Busco outliers y absurdos y los paso a missing ----
# En el summary vimos que había una edad de 144
datos$age<- replace(datos$age, which(datos$age > 100), NA)
datos$loan_amnt[abs(scores(datos$loan_amnt, type="z"))>4] # Aparece el valor 35000 repetido unas 160 veces
##   [1] 35000 35000 35000 35000 35000 35000 35000 35000 35000 35000 35000
##  [12] 35000 35000 35000 35000 35000 35000 35000 35000 35000 35000 35000
##  [23] 35000 35000 35000 35000 35000 35000 35000 35000 35000 35000 35000
##  [34] 35000 35000 35000 35000 35000 35000 35000 35000 35000 35000 35000
##  [45] 35000 35000 35000 35000 35000 35000 35000 35000 35000 35000 35000
##  [56] 35000 35000 35000 35000 35000 35000 35000 35000 35000 35000 35000
##  [67] 35000 35000 35000 35000 35000 35000 35000 35000 35000 35000 35000
##  [78] 35000 35000 35000 35000 35000 35000 35000 35000 35000 35000 35000
##  [89] 35000 35000 35000 35000 35000 35000 35000 35000 35000 35000 35000
## [100] 35000 35000 35000 35000 35000 35000 35000 35000 35000 35000 35000
## [111] 35000 35000 35000 35000 35000 35000 35000 35000 35000 35000 35000
## [122] 35000 35000 35000 35000 35000 35000 35000 35000 35000 35000 35000
## [133] 35000 35000 35000 35000 35000 35000 35000 35000 35000 35000 35000
## [144] 35000 35000 35000 35000 35000 35000 35000 35000 35000 35000 35000
## [155] 35000 35000 35000 35000 35000 35000 35000 35000
datos$annual_inc[abs(scores(datos$annual_inc, type="z"))>4] # El tercer cuartil está a mucha distancia del máximo
##   [1]  365000  350000  408000  400000  480000  500000  357000  385000
##   [9]  400000  350000  395000 1782000  900000  750000  350000  400000
##  [17]  350000  480000  336000  372000  325000  381450  370000  350000
##  [25]  487000  648000  390000  350000  410000  450000  336000  552000
##  [33]  430000  334000  480000  330000  500000  780000  500000  948000
##  [41] 1200000  408000  350000  445000  828000  330000 1900000  360000
##  [49]  550000  667680  564000  889000  550000 1200000  900000  400000
##  [57]  528000  351700  741600  840000 6000000  660000  415000  648000
##  [65]  612000  334000  504000  500000  425000  800000  340000  600000
##  [73]  510000  420000  325000  384000  600000  636000 1200000  780000
##  [81]  444000  762000  708000  371000 1440000  900000  900000 1362000
##  [89]  348000  720000  480000  360000  780000  700000  322400  450000
##  [97]  600000  700000  390000  350000  500000 2039784  572400  375000
## [105]  522000  616000  480000
# Miro las distribuciones que parecen sospechosas
ggplot(datos, aes(x=annual_inc)) +
  geom_density(alpha=.3, show.legend = FALSE) +
  theme(axis.title.x=element_blank()) + ggtitle("Dist. Ingresos anuales")

# Ambas son aceptables

## 1.3 Sustituyo los datos missing ----
table(is.na(datos$impago)) # No tiene
## 
## FALSE 
## 29092
table(is.na(datos$loan_amnt)) #No tiene
## 
## FALSE 
## 29092
table(is.na(datos$int_rate)) # Tiene, los sustituyo por el valor de préstamos similares
## 
## FALSE  TRUE 
## 26316  2776
datos<- rellenaMissings(datos, "int_rate", "loan_amnt", FALSE)
table(is.na(datos$Calidad_Crediticia))
## 
## FALSE 
## 29092
table(is.na(datos$age)) # Se que solo hay un missing en edad, lo sustrituyo por la media
## 
## FALSE  TRUE 
## 29091     1
datos[is.na(datos$age),"age"]<- ceiling(mean(datos$age, na.rm=T)) 
table(is.na(datos$emp_length)) # Tiene, los sustituyo por valores de personas similares
## 
## FALSE  TRUE 
## 28283   809
datos<- rellenaMissings(datos, "emp_length", "age", TRUE)
table(is.na(datos$home_ownership)) # No tiene
## 
## FALSE 
## 29092
table(is.na(datos$annual_inc)) # No tiene
## 
## FALSE 
## 29092

# 2. Tratamiento  variables ----
summary(datriesgos)
##      impago         loan_amnt        int_rate     Calidad_Crediticia
##  Min.   :0.0000   Min.   :  500   Min.   : 5.42   A:9649            
##  1st Qu.:0.0000   1st Qu.: 5000   1st Qu.: 7.90   B:9329            
##  Median :0.0000   Median : 8000   Median :10.99   C:5748            
##  Mean   :0.1109   Mean   : 9594   Mean   :11.00   D:3231            
##  3rd Qu.:0.0000   3rd Qu.:12250   3rd Qu.:13.47   E: 868            
##  Max.   :1.0000   Max.   :35000   Max.   :23.22   F: 211            
##                                   NA's   :2776    G:  56            
##    emp_length      home_ownership    annual_inc           age       
##  Min.   : 0.000   MORTGAGE:12002   Min.   :   4000   Min.   : 20.0  
##  1st Qu.: 2.000   OTHER   :   97   1st Qu.:  40000   1st Qu.: 23.0  
##  Median : 4.000   OWN     : 2301   Median :  56424   Median : 26.0  
##  Mean   : 6.145   RENT    :14692   Mean   :  67169   Mean   : 27.7  
##  3rd Qu.: 8.000                    3rd Qu.:  80000   3rd Qu.: 30.0  
##  Max.   :62.000                    Max.   :6000000   Max.   :144.0  
##  NA's   :809
summary(datos) # La inputación no ha modificado significativamente los datos
##  impago      loan_amnt        int_rate     Calidad_Crediticia
##  0:25865   Min.   :  500   Min.   : 5.42   A:9649            
##  1: 3227   1st Qu.: 5000   1st Qu.: 8.00   B:9329            
##            Median : 8000   Median :10.99   C:5748            
##            Mean   : 9594   Mean   :11.01   D:3231            
##            3rd Qu.:12250   3rd Qu.:13.24   E: 868            
##            Max.   :35000   Max.   :23.22   F: 211            
##                                            G:  56            
##    emp_length      home_ownership    annual_inc           age      
##  Min.   : 0.000   MORTGAGE:12002   Min.   :   4000   Min.   :20.0  
##  1st Qu.: 2.000   OTHER   :   97   1st Qu.:  40000   1st Qu.:23.0  
##  Median : 4.000   OWN     : 2301   Median :  56424   Median :26.0  
##  Mean   : 6.242   RENT    :14692   Mean   :  67169   Mean   :27.7  
##  3rd Qu.: 8.626                    3rd Qu.:  80000   3rd Qu.:30.0  
##  Max.   :62.000                    Max.   :6000000   Max.   :94.0  
## 
varObjBin<- datos$impago

En primer lugar debo señalar que el dataset tenía muy pocos NAs, para inputarlos lo que he hecho es rellenar con valores que tienen otra variable que yo he considerado importante, similar, añadiendo un término aleatorio, así el tipo de interés según datos con cantidades prestadas similares, y la duración del empleo con datos similares en edad. Para la edad como solo faltaba un dato lo he reemplazado con la media.

El siguiente paso va a consisstir en hacer variables categóricas de todas las contínuas. Voy a representar la V de Cramer de las variables frente a impago para seleccionar dos grupos de las más significativas (que usaré más tarde para hacer modelos) y voy ha hacer algunos gráficos representativos de la distribución de las variables numéricas más significativas, para comprender algo mejor los datos. Por último voy a hacer algunas transformaciones sobre las variables cuantitativas y ha separar los datos en muestra de entrenamiento y muestra de test.

## 
## $int_rate

## 
## $Calidad_Crediticia

## 
## $emp_length

## 
## $home_ownership

## 
## $annual_inc

## 
## $age

# breaks_adj<- woebin_adj(datos, age, bins)
# loan_amnt=c("3000", "5500", "7000", "11000", "14500", "16500", "21500"), 
# int_rate=c(6.5, 8, 11, 15), 
# Calidad_Crediticia=c("A", "B", "C", "D%,%E%,%F%,%G"), 
# emp_length=c(1.5, 2.5, 6.5, 10.5, 13), 
# home_ownership=c("MORTGAGE", "OTHER%,%OWN", "RENT"), 
# annual_inc=c(40000, 65000, 95000), 
# age=c(23, 24, 38)
# Incorporo las variables binneadas y agrupadas
datos$loan_amnt.fac<- cut(datos$loan_amnt, breaks=c(0, 3000, 5500, 7000, 11000, 14500, 16500, 21500, Inf))
datos$int_rate.fac<- cut(datos$int_rate, breaks=c(0, 6.5, 8, 11, 15, Inf))
datos$emp_length.fac<- cut(datos$int_rate, breaks=c(0, 6.5, 8, 11, 15, Inf))
datos$annual_inc.fac<- cut(datos$annual_inc, c(0, 40000, 65000, 95000, Inf))
datos$age.fac<- cut(datos$age, c(0, 23, 24, 38, Inf))
datos$home_ownership.gro<- droplevels(replace(datos$home_ownership, which(datos$home_ownership == "OWN"), "OTHER"))
datos$Calidad_Crediticia.gro<- droplevels(replace(datos$Calidad_Crediticia,
                                                  which(datos$Calidad_Crediticia %in% c("E", "F", "G")), "D"))
## 2.3 Importancia de las variables ----
# Calcula el V de Cramer
CrossTable(varObjBin) # 88.9% de ceros (NO impago), 11.1% de unos (SI impagos)
## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  29092 
## 
##  
##           |         0 |         1 | 
##           |-----------|-----------|
##           |     25865 |      3227 | 
##           |     0.889 |     0.111 | 
##           |-----------|-----------|
## 
## 
## 
## 
par(mar = c(4.1,10,4.1,4.1))
# Dejo el impago para asegurarme de la calidad de la función
graficoVcramer(datos, datos$impago, "Influencia de las variables")

Mirando los gráficos si encontramos que el tipo de interés depende de la duración del empleo y de la Calidad Crediticia, mientras que la cantidad prestada muestra mayor variación frente a los ingresos anuales. Las variables

Debo de decir que no termino de entender el uso de Welch test aquí ya que según he leído se emplea para constrastar similitud entre distribuciones normales, y nosotros estamos empleando una de las variables categoóricas. No obstante no he tenido el tiempo que me hubiera gustado para informarme sobre el test. Así que en la variable para la que el p-valor es dudoso, he probado una transformación y como lo mejora la he incorporado a los datos.

El siguiente paso va a consistir ya en diseñar los modelos, voy a hacer 7 uno con todas las variables, dos con las listas de variables seleccionadas, dos con el criterio de selección de variables AIC, primero en las dos direcciones y luego sólo hacia atrás, y otro modelo haciendo esto mismo con el criterio de selección BIC.

De momento los modelos ya están entrenados. Vemos que pese a las dummies el número de variables no se ha disparado, lo interesante es ahora hacer validación cruzada repetida y ver cuál es el modelo que mejor se comporta

Lo primero que conviene notar es que todos los modelos presentan una área bajo la curva ROC de muy similar, entre el menor y el mayor valor hay tan solo 0.0045189. A priori lo mejor sería coger el de menor número de variables, sin embargo, yo he seleccionado el modelo mAIC.b que con 21 variables mejora en 4 milésimas a aquellos que presentan 10 variables menos y a todos los que presentan más variables. Tal vez 10 variables para 4 milésimas no merecen la pena, pero 21 variables no son demasiadas.

Podemos ver que el PSI es: 3.610^{-4}, suficientemente bajo como para dar por exitoso el modelo de scoring, que presenta además una prescisión sobre el set de test de 86.855% que es bastante elevada, así que debemos darnos por satisfechos.

Muchas gracias

Pese a que el set de ejercicios me ha parecido desproporcionadamente largo, creo que el esfuerzo ha merecido la pena, gracias por su trabajo.