Construcción de curvas suaves con series de Fourier

  1. Seleccionar cinco días consecutivos completos (que tengan las 24 mediciones en cada día) de la base de datos de \(NO_2\).

  2. Usar las primeras siete funciones de la \(CONS\) de cosenos para ajustar un estimador de series en cada uno de los cinco días elegidos, usando la hora como variable \(x\) (convertir las horas a valores en \([0,1]\)) y los niveles horarios de \(NO_2\) como respuesta.

  3. Representar las cinco suavizaciones en un solo gráfico y comentar los resultados.

Librerias necesarias

library(readxl)
library(dplyr)
library(lubridate)
library(ggplot2)
library(purrr)

Carga de Datos

Para efectos de este ejercicio Se hace uso del dataset \(DatosLaFlora\) el cual contiene información sobre mediciones horarias de \(NO_2\) y \(O_3\) en (\(\mu g/m3\)) y sobre algunas variables de clima en el norte de Cali(Colombia). Se dispone además de un campo que informa sobre el día y la hora a la que se tomó la medición.

Datos <- read_excel(ruta)

Ajuste de los datos

Lo que se busca en esta sección es conformar un dataframe con las variables (\(NO_2\), \(dia\) y \(hora\_normada\)) del dataset \(DatosLaFlora\), las variables \(dia\) y \(hora\_normada\) se obtienen a partir de la variable \(Fecha \& Hora\); \(hora\_normada\) se usa para restringimos al intervalo \([0,1]\), este espacio de las funciones cuadrado integrables se denota en general como \(L_2[0,1]\) el cual satisface las propiedades generales de un espacio vectorial.

La transformación se realiza de la siguiente forma:

para cada hora \(i \in [1,24]\) se aplica la ecuación para normalizar:

\[x_i=\frac{(2_i-1)}{2n}, \hspace{1 cm} \forall i=1,2,...,24\]

x <- Datos  %>% 
  dplyr::select(`Fecha & Hora`, NO2) %>% 
  mutate(dia = day(`Fecha & Hora`)) %>%
  filter(dia %in% c(9,10,11,12,13)) %>%  
  group_by(dia) %>%
  
  mutate(hora = row_number()) %>% 
  ungroup() %>%
        
  mutate(hora_normada = (2*hora - 1)/(2*24)) %>%
  mutate(dia = as.factor(dia))

Obtenemos entonces el siguiente \(tibble\):

## # A tibble: 120 x 5
##    `Fecha & Hora`        NO2 dia    hora hora_normada
##    <dttm>              <dbl> <fct> <int>        <dbl>
##  1 2011-08-09 00:00:00  66.9 9         1       0.0208
##  2 2011-08-09 01:00:00  57.9 9         2       0.0625
##  3 2011-08-09 02:00:00  54.5 9         3       0.104 
##  4 2011-08-09 03:00:00  47.2 9         4       0.146 
##  5 2011-08-09 04:00:00  39.8 9         5       0.188 
##  6 2011-08-09 05:00:00  36.1 9         6       0.229 
##  7 2011-08-09 06:00:00  36.7 9         7       0.271 
##  8 2011-08-09 07:00:00  37.1 9         8       0.312 
##  9 2011-08-09 08:00:00  37.6 9         9       0.354 
## 10 2011-08-09 09:00:00  41.7 9        10       0.396 
## # ... with 110 more rows

Exploración de los datos

En esta sección se genera un solo gráfico de línea con las mediciones de \(NO_2\) vs \(hora\_normada\) para cada uno de los días seleccionados.

ggplot() +
  geom_line(data = x, aes(x = hora_normada, y = NO2, color = dia))+
  labs(y = expression(NO[2]))

Y adicionalmente gráficos de puntos individuales para cada día.

ggplot()+
  geom_point(data = x, aes(x = hora_normada, y = NO2)) +
  labs(y = expression(NO[2]))+ facet_wrap(~dia)

Podemos observar y comparar el comportamiento de los datos a lo largo de los días tanto de forma conjunta como individual.

Ajuste del modelo

El propósito de este trabajo es ajustar un modelo básico al conjunto de observaciones de 5 dias, usando el siguiente estimador de series de Furier:

\[\mu_{\lambda}(x)=\sum_{j=1}^\lambda \beta_{\lambda j}f_j(x)\]

Haciendo uso de una sucesión ortonormal completa \((CONS)\) de funciones, la cual es la siguiente base de cosenos:

\[\sqrt{2}cos[(j-1)\pi x], \hspace{1 cm} \forall j=2,\dots \]

Se ajusta entonces un modelo mediante transformaciones de Furier y se toma arbitrariamente un valor\(\lambda = 7\).

\[\mu_{\lambda}(x)=\sum_{j=1}^\lambda \beta_{\lambda j}\sqrt{2}cos[(j-1)\pi x]\]

Construcción data frame

En esta sección deseamos construir un dataframe tomando como variable respuesta \(NO_2\) y como variables predictoras \(f_1, f_2, \dots, f_j\) donde \(j=\lambda\) y \(f\) es la base de cosenos \((CONS)\) que elegimos previamente.

base_cons <- function(x,j){
              sqrt(2)*cos((j-1)*pi*x)
              }

df <-  x %>% 
  dplyr::select(NO2) 

hora_normada <- x$hora_normada; lambda <- 7; i <- NULL; aux <- NULL; y <- NULL

for (i in 2:lambda) {
  aux <- base_cons(hora_normada, i)
  y[[i]] <- aux
  }

y <- data.frame(matrix(unlist(y), ncol = lambda-1)) # Variables predictoras
df <-  bind_cols(df, y)

Una vez construido el data frame procedemos con la estimación por mínimos cuadrados de los coeficientes de Furier haciendo uso de la función lm(). Esta estimación debe realizarse para cada uno de los 5 días.

lm_grupo <- function(x){
  model <- lm(NO2 ~  X1+X2+X3+X4+X5+X6, data = x)
  fitted_NO2 <- predict(model) %>% 
    as_tibble
  bind_cols(x, fitted_NO2)
}

df <- df  %>% 
  mutate(Dia = as.factor(x$dia)) %>%
  group_split(Dia) %>% 
  map(.x = ., .f = lm_grupo) %>% 
  bind_rows()

Obtenemos entonces los fitted values del modelo para cada día y los agregamos al final del dataframe df.

fitted <- as.data.frame(df$value)
names (fitted) = "fitted"
df = bind_cols(x, fitted)
## # A tibble: 120 x 6
##    `Fecha & Hora`        NO2 dia    hora hora_normada fitted
##    <dttm>              <dbl> <fct> <int>        <dbl>  <dbl>
##  1 2011-08-09 00:00:00  66.9 9         1       0.0208   64.1
##  2 2011-08-09 01:00:00  57.9 9         2       0.0625   60.7
##  3 2011-08-09 02:00:00  54.5 9         3       0.104    54.7
##  4 2011-08-09 03:00:00  47.2 9         4       0.146    47.5
##  5 2011-08-09 04:00:00  39.8 9         5       0.188    40.7
##  6 2011-08-09 05:00:00  36.1 9         6       0.229    35.9
##  7 2011-08-09 06:00:00  36.7 9         7       0.271    34.0
##  8 2011-08-09 07:00:00  37.1 9         8       0.312    35.5
##  9 2011-08-09 08:00:00  37.6 9         9       0.354    40.0
## 10 2011-08-09 09:00:00  41.7 9        10       0.396    46.4
## # ... with 110 more rows

Visualización

Finalmente con el nuevo dataframe df podemos generar un gráfico donde se puedan observar los datos reales (puntos) y los datos ajustados por el modelo (línea) de la variable \(NO_2\) vs \(hora\_normada\) para cada uno de los días.

ggplot()+
  geom_point(data = df, aes(x = hora_normada, y = NO2)) +
  geom_line(data = df, aes(x =hora_normada, y = fitted)) +
  labs(y = expression(NO[2])) +
  labs(title = "Estimacion con base de cosenos") +
  labs(subtitle = expression(lambda==7)) +
  facet_wrap(~dia)

Tenemos entonces la estimación del comportamiento de \(NO_2\) para los días \(9,10,11,12\) y \(13\) usando series de Furier con base de cosenos y con un \(\lambda=7\), podríamos decir que \(\hat{\mu}_7\) es una buena aproximación a \(\mu\), pero dado que \(\lambda\) fue seleccionado de forma arbitraria no podemos asegurar que \(\lambda=7\) es la mejor opción, ¿cuál valor de \(\lambda\) seria entonces “la mejor” elección? de esto precisamente nos ocuparemos en la siguiente entrega.

El ajuste de este modelo se llevo acabo basándose en material expuesto en (Olaya, 2012).

Bibliografía

  • Olaya, J. (2012). Métodos de Regresión No Paramétrica. Universidad del Valle.

  • R Core Team. (2013). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing. http://www.r-project.org/