May 10, 2017

Datos

Fuente: http://lib.stat.cmu.edu/datasets/fraser-river. - \(946\) observaciones de una variable (nivel medio mensual).

river <- read.table('fraser-river.txt'); head(river)
##     V1
## 1  485
## 2 1150
## 3 4990
## 4 6130
## 5 4780
## 6 3960

Datos

Objetivo

Parte \(A\)

  • Se pide estimar la probabilidad de que el río desborde en un mes superando los \(10000 m^3/s\).

Modelo

Parte \(A\cdot B\)

  • Modelamos el nivel del río utilizando un PAR(1) - Periodic Autoregressive Model of order 1.
  • Un modelo AR(1) tiene la forma

\[ y_t \sim Normal(\phi_0+\phi_1 y_{t-1}, \varphi)\]

  • En un modelo PAR(1) incluimos componentes que varían según la estación. En nuestro caso el período de la serie (el número de meses) es 12.

\[ y_t \sim Normal(\phi_{0s}+\phi_{1s} y_{t-1}, \varphi_s), s=1,2, \ldots 12.\]

Modelo Jerárquico

  • \(Y_{t} \big|\{ \phi_{0s}, \phi_{1s}, \varphi_s, Y_{t-1} \} \sim Normal(\phi_{0s} + \phi_{1s} Y_{t-1}, \varphi_s)\),
  • \(\phi_{0s}\sim Normal(\mu_0, \sigma_0^2)\),
  • \(\phi_{1s}\sim Normal(\mu_1, \sigma_1^2)\),
  • \(\varphi_{s}^2\sim InvGamma(\alpha, \beta)\).

Estamos interesados en \(P(Y_t > 10000)\). Note que \(t\) depende de \(s\), por ejemplo,

\[ t = 12 \cdot n + s, \quad n \in \mathbb{Z}_+.\]

Datos transformados y centrados

Estos son nuestros datos. Aplicamos una transformación logarítmica.

Implementación en stan

Código completo

model{
  // Formulación del modelo
  int season;
  season = 2;
  for(i in 2:N){
    y[i] ~ normal( phi0[season]+phi1[season]*y[i-1], varphi[season] );
    season=season+1;
    if(season>12){
      season=1;
    }
  }
  
  for (s in 1:12){
    phi0[s] ~ normal( mu0, sigma0 ); // densidad a-priori normal
    phi1[s] ~ normal( mu1, sigma1 );
    varphi2[s] ~ inv_gamma( alpha, beta ); // densidad a-priori inversa gamma
  }
}

Implementación en stan

generated quantities{
  // Valores simulados para estimar

  vector[N_new] y_new;
  int season;
  season = 2;

  y_new[1]=jan_new;

  for(i in 2:N_new){
    y_new[i] = normal_rng( phi0[season]+phi1[season]*y_new[i-1], varphi[season] );
    season=season+1;
    if(season>12){
      season=1;
    } 
  }
}

Trayectorias Simuladas

Densidad de probabilidad marginal

Probabilidad de desbordamiento por mes

April May June July August September
Prob. 0 0.0017333 0.0472 0.0145333 0.0002667 0

Gráficos de diagnóstico

Diagnóstico usando el traceplot del muestreo de Gibbs para la variable \(\phi_{01}\).

Resultados

Parte \(B^{-1}\cdot A \cdot B\).

  • La probabildad estimada de que el río desborde superando los 10000 \(m^3/s\) en un mes es de \(0.0053\).
  • La probabilidad de que el río desborde en el mes de Junio es de \(0.0472\).
mean(result.sim$y_new[,]>log(10000))  
mean(result.sim$y_new[,12*(2:11)+6]>log(10000))

Muchas Gracias!