Repositorio Institucional de Documentos

Trabajo de Fin de Grado Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín Autora Carolina Pelegrín Cuartero Directores Jesús Asín Lafuente Jesús Fernández Cascán Ponente Beatriz Lacruz Casaucau Ciencias Ambientales Escuela Politécnica Superior de Huesca Junio de 2014 ÍNDICE 1. INTRODUCCIÓN ............................................................................................................................. 1 2. DATOS DISPONIBLES Y ANÁLISIS EXPLORATORIOS ........................................................................ 7 2.1 Análisis y pretratamiento de la base de datos ............................................................................ 7 2.1.1 Análisis de las bases de datos disponibles .............................................................................. 7 2.1.2 Pretratamiento de la bases de datos observados .................................................................. 8 2.2 Estudio de los regímenes de vientos, real y predicho ................................................................ 10 2.3 Análisis de correlación .............................................................................................................. 11 3. METODOLOGÍA. CONSTRUCCIÓN DE MODELOS DE REGRESIÓN .................................................. 12 3.1 Procedimiento de construcción de modelos candidatos ........................................................... 13 3.1.1 Estudio de variables potencialmente útiles. Construcción de la base de datos ................... 13 3.1.2 Estudio de una posible transformación de variables ............................................................ 15 3.1.3 Estrategia de construcción de modelos ................................................................................ 15 3.1.4 Criterios de bondad de ajuste y selección de variables ........................................................ 16 3.2 Selección de modelos candidatos .............................................................................................. 17 3.3 Estudio de la adecuación de los modelos seleccionados ........................................................... 18 3.4 Estudio de la capacidad predictiva del modelo ......................................................................... 19 3.5 Estudio del funcionamiento operativo de los modelos seleccionados. Realización de pronósticos ............................................................................................................................................. 20 3.6 Selección del modelo óptimo y estudio operativo ..................................................................... 22 4. RESULTADOS ............................................................................................................................... 23 4.1 Análisis exploratorio.................................................................................................................. 23 4.1.1 Análisis, pretratamiento y selección de la base de datos ..................................................... 23 4.1.2 Estudio del régimen de vientos en Bailín y de la adecuación del modelo de predicciones de AEMET …………………………………………………………………………………………………………………………………………25 4.1.3 Análisis de correlación .......................................................................................................... 33 4.1.4 Principales conclusiones del análisis exploratorio ................................................................ 36 4.2 Desarrollo de modelos de regresión .......................................................................................... 36 4.2.1 Estimación del modelo de regresión lineal simple para los datos de predicción existentes 37 4.2.2 Desarrollo de modelos de regresión candidatos a ser el óptimo buscado ........................... 41 4.2.3 Selección de modelos candidatos ......................................................................................... 45 4.2.4 Adecuación de los modelos seleccionados ........................................................................... 50 4.2.5 Estudio del funcionamiento operativo de los modelos seleccionados. Pronóstico ............. 54 4.3 Selección del modelo óptimo, estudio de su capacidad operativa e implementación del procedimiento ........................................................................................................................................ 64 5. CONCLUSIONES ........................................................................................................................... 66 Agradecimientos .................................................................................................................................. 69 Bibliografía .......................................................................................................................................... 70 Anexo I. Figuras ................................................................................................................................... 73 Análisis exploratorio ............................................................................................................................... 73 Desarrollo de modelos de regresión ....................................................................................................... 74 Anexo II. Conceptos estadísticos de interés ......................................................................................... 99 Anexo III. Estudios adicionales ........................................................................................................... 107 RESUMEN Este Trabajo Fin de Grado se desarrolla en el marco del proceso de desmantelamiento del vertedero de residuos organoclorados de Bailín (Huesca) y su transferencia a un depósito seguro, llevado a cabo durante el verano de 2014. El fin último del trabajo ha sido crear una herramienta que permita anticipar posibles eventos dispersivos de la contaminación provocados por el viento, capaces de agravar la problemática ambiental existente en la zona. Para ello, en primer lugar se ha realizado un análisis exploratorio del régimen de vientos en el vertedero de Bailín para el periodo de verano de 2012, así como un estudio que demuestra que el modelo de predicciones de AEMET existente no resulta útil de forma directa. Por ello, se han desarrollado modelos estadísticos de predicción de la velocidad del viento a 24 horas en forma de modelos de regresión (MR) que han permitido el downscaling del modelo de predicciones de AEMET, adaptándolo a las condiciones específicas del área y el periodo de estudio. El proceso de construcción de MR se ha llevado a cabo en el entorno del software estadístico R mediante los siguientes pasos: 1) Estudio de parámetros potencialmente útiles (de predicción AEMET y otros), y construcción de una base de datos con ellos; 2) Desarrollo de una estrategia de construcción y selección de modelos (sin y con interacciones), haciendo uso de métodos stepwise de selección de parámetros, así como de criterios de selección de modelos como bondad de ajuste, desviación típica residual de validación cruzada o realización de test ANOVA; 3) Estudio de la adecuación de los modelos mediante el análisis de residuos; 4) Estudio del funcionamiento de los modelos candidatos finales; 5) Selección del modelo óptimo e implementación para su uso operativo. El modelo finalmente seleccionado ha aumentado el porcentaje de explicación de la variabilidad de la velocidad del viento hasta un 57%, respecto al 14% obtenido con las predicciones de AEMET. La metodología creada también pretende servir de base para la resolución de problemas de predicción similares. ABSTRACT This study was developed within the framework of the Bailín landfill (Huesca, Spain) dismantlement and its waste transfer to a new and safe location. The aim of the work has been to create a tool that will allow to anticipate wind events, which could produce pollution dispersion, aggravating the environmental situation at the area. First, an exploratory analysis of the wind regime at the area has been done. This first studies showed that the AEMET (Agencia Estatal de Meteorología) prediction model didn’t reflect the real wind speed at the area and period of interest (Bailín, summer of 2012). That’s the reason why a regression model (RM) has been developed; it has allowed us to downscale the AEMET wind speed predictions, adapting them to the real wind conditions in Bailín. The model provides the wind speed prediction at the Bailín area, 24 hours in advance. The construction of the RM was carried out within the R environment, through the following steps: 1) Study of the useful parameters to be included at the model, and construction of a database with them; 2) Developing of a model-construction and model-selection strategy (with and without interactions), using stepwise methods, as well as model selection criteria (R 2 adjusted, residual deviation of cross validation or ANOVA test); 3) Study of the models’ fit, by studying its residuals; 4) Study of the finally-selected-models’ running at the area and period of interest; 5) Final best-model selection. The model finally selected has been able to increase the explanation of the real wind speed variability up to a 57%, with regard to the 14% explained by the AEMET predictions. The methodology developed also pretends to provide a tool for solving similar prediction problems. 1 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín 1. INTRODUCCIÓN HCH TÉCNICO Y LINDANO Desde la antigüedad, la práctica de la agricultura ha estado unida al uso de compuestos destinados a mejorar la producción y evitar las plagas. A comienzos del siglo XX, y especialmente tras la Segunda Guerra Mundial, comenzó la utilización masiva de productos químicos organoclorados 1 como pesticidas, tanto en la agricultura como en el control de parásitos (Tadeo, 2008). Uno de los que se introdujo más rápidamente en el mercado mundial de productos fitosanitarios fue el hexaclorociclohexano (HCH), gracias a sus propiedades insecticidas, su rápido efecto y su fácil obtención industrial (Li et al., 1998). El HCH es un insecticida organoclorado de amplio espectro que se usó principalmente desde la Segunda Guerra Mundial hasta la década de 1990 (Forter, 1995; Breivik et al., 1999). Su nomenclatura es 1,2,3,4,5,6-hexaclorociclohexano, pues todos sus hidrógenos se encuentran sustituidos por átomos de cloro (ver Figura 1 izda.). Lo componen 8 isómeros diferentes, y durante los primeros años, se utilizó el conjunto de todos ellos (conocido como HCH técnico) como un insecticida económico. Más adelante, se descubrió que su actividad insecticida radicaba únicamente en el isómero gamma (γ-HCH, ver Figura 1 dcha.), el cual comenzó a refinarse a partir del HCH técnico, y se comercializó con el nombre de lindano (Walker et al., 1999). El lindano representa únicamente entre el 8% y el 15% del HCH técnico (Vijgen et al., 2011), por lo cual su fabricación genera una gran cantidad de residuos, formados por el resto de isómeros “inactivos”, principalmente α-HCH, β-HCH y δ-HCH (Bodenstein, 1972). Figura 1. Composición química del Hexaclorociclohexano (HCH técnico, izda.) y del lindano (γ-HCH, dcha.), el isómero que posee las propiedades insecticidas. (Fernández, 2004). Los compuestos derivados del HCH pertenecen al grupo de los pesticidas persistentes por su baja biodegradabilidad, la cual, junto a su toxicidad, supone una grave amenaza tanto para el medio ambiente como para la salud (Willett et al., 1998). Sin embargo, durante muchos años se 1 Compuestos químicos orgánicos, es decir, con un esqueleto de átomos de carbono, en el cual alguno de los átomos de hidrógeno han sido sustituidos por átomos de cloro. 2 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín consideró que dichos compuestos eran insolubles en el agua e inocuos en general, y fueron vertidos al medio ambiente sin ningún control e incluso utilizados en la construcción (Forter, 1995; Vijgen 2006) o pavimentación (Torres et al., 2012). Esto provocó importantes focos de contaminación a nivel mundial, los cuales se han detectado a lo largo de los últimos años (Götz et al., 2013; Jit et al., 2010; Vijgen et al., 2011; Weber and Varbelow 2012; Wycisk et al., 2012). En la década de 1970 comenzaron a detectarse ciertos problemas ambientales ocasionados por los derivados del HCH, lo que propició su prohibición en muchos países (Vijgen et al., 2011; Voldner and Li, 1995; Walker et al., 1999) Sin embargo, no fue hasta el año 2009, cuando los isómeros α-HCH, β-HCH y γ-HCH (lindano) se incluyeron en la Convención de Estocolmo sobre Contaminantes Persistentes Orgánicos (COPs), cuyo principal objetivo es el de proteger la salud humana y el medio ambiente frente a este tipo de contaminantes. Desde entonces, el problema de la contaminación derivada del uso y la producción del HCH y el lindano, junto con sus depósitos de residuos, se reconoció como un asunto de preocupación a nivel mundial (Vijgen et al., 2011) debiendo identificarse los lugares contaminados, así como adoptar las medidas necesarias para gestionar dichos residuos. En resumen, la producción de lindano ha generado el mayor almacén mundial de COPs, estimado entre 4.8 y 7.2 millones de tm (Vijgen, 2006; Vijgen et al., 2011). El tratamiento de las áreas contaminadas debe ser abordado como parte de la implementación de la Convención de Estocolmo. ÁREA DE ESTUDIO. EL BARRANCO DE BAILÍN En España, uno de los principales focos de contaminación por lindano se localiza en los alrededores de Sabiñánigo (Huesca), donde se calcula que la empresa INQUINOSA (Industrias Químicas del Noroeste, ver Figura 2), que operó en la zona entre 1975 y 1988, vertió en los alrededores 115 000 tm de residuos peligrosos organoclorados, en gran parte isómeros de HCH (Vijgen, 2006; Fernández et al., 2012) sin ningún tipo de política medioambiental. Uno de los principales puntos de depósito fue el barranco de Bailín, nuestra zona de estudio (ver Figura 3). Según Fernández et al. (2012), en él se vertieron, entre los años 1984 y 1988, entre 30 000 y 80 000 tm 3 de residuos sólidos de HCH, así como 2000 toneladas de residuo líquido DNAPL (dense nonaquous phase liquid, ver Figura 4). Estas estimaciones lo sitúan como uno de los mayores depósitos de HCH a nivel mundial (Vijgen, 2006). En el año 2005 se descubrió que el DNAPL estaba percolando hacia el río Gállego, que se encuentra a tan solo 800 metros. Por ello, en 2007 se decidió el desmantelamiento y transferencia del vertedero de Bailín a otro emplazamiento que tuviese las medidas de seguridad oportunas. Como nuevo emplazamiento, y teniendo en cuenta preocupaciones locales y consideraciones de seguridad, se eligió un sitio cercano en el mismo valle, el cual se ha construido con importantes medidas de aislamiento (ver Figura 5). La capacidad del nuevo vaso es de 250 000 m 3 , suficiente para almacenar el contenido del vertedero de Bailín, los residuos que hay dentro de INQUINOSA y los suelos afectados alrededor del vertedero de Sardás, el segundo punto de mayor depósito de residuos, después de Bailín (Fernández et al., 2012). 3 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín Figura 2. Antigua fábrica de INQUINOSA, hoy en día abandonada. Fotografía extraída del Diario del Alto Aragón. Figura 3 Localización geográfica de los puntos de contaminación por residuos de HCH alrededor de Sabiñánigo, incluido el vertedero de Bailín (Bailín dump) (Fernández et al., 2012). Figura 4. Residuo líquido DNAPL en el área del vertedero de Bailín. Fotografía de Alfonso Pardo. Figura 5. Nuevo vaso construido para albergar los residuos organoclorados del antiguo vertedero de Bailín. Fotografía de Jesús Fernández. Según Fernández et al. (2012), la construcción de dicha infraestructura comenzó en 2010 y las obras de desmantelamiento y transferencia se prevén para el verano de 2014, entre los meses de mayo a octubre, cuando las precipitaciones son menores. Aunque esto no supone una eliminación final de los residuos de HCH, el desmantelamiento es necesario dado el impacto negativo que las percolaciones de DNAPL suponen para el río Gállego. Según lo dispuesto en la Convención de Estocolmo, en un futuro se deberán eliminar los residuos almacenados en el nuevo vaso. En Bailín, este proceso se llevará a cabo cuando esté disponible un tratamiento asequible. 4 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín DISPERSIÓN DE PARTÍCULAS EN VERTEDEROS Y PROTOCOLO DE ACTUACIÓN EN BAILÍN Hasta hace no muchos años no existían normas sobre deposición de residuos en vertederos y, tanto empresas como núcleos urbanos, depositaban sus residuos en los alrededores de su actividad sin ningún tipo de control. Como ya se ha mencionado anteriormente, lo mismo ocurrió con los residuos de fabricación del lindano. Hoy en día, el principal método de eliminación de residuos en Europa sigue siendo su deposición en vertederos (Giusti, 2009), pero la UE controla esta actividad gracias a reglas muy estrictas en su legislación (Directiva 2008/98/EC). Se espera que durante el proceso de desmantelamiento y transferencia de los residuos contenidos en el vertedero de Bailín al nuevo vaso, se emita una gran variedad de contaminantes al aire, principalmente gases y partículas (Chalvatzaki et al., 2010; Kopanakis et al., 2011), a causa de los procesos que se llevarán a cabo de descomposición, resuspensión y operaciones típicas de vertedero, como carga y descarga de residuos, o transporte de los mismos. Estas emisiones estarán también influenciadas por las condiciones climáticas locales, especialmente por la presencia e intensidad del viento (Chalvatzaki et al., 2012). Por ello, las operaciones llevadas a cabo en el vertedero de Bailín se suspenderán cuando se prevean lluvias e intensidades de viento elevadas (Fernández et al., 2012). En resumen, uno de los mayores retos en el desmantelamiento y transferencia de los residuos contenidos en el vertedero de Bailín, consiste en evitar la dispersión atmosférica de partículas de lindano en la medida de lo posible, por los riesgos ambientales y para la salud que ésta supondría (Pope et al., 1995). Para ello, y según el protocolo de actuación de las obras, éstas se detendrán cuando se prevean vientos superiores a 40 km/h. Esta parada no será instantánea, es decir, una vez que esté sucediendo o haya sucedido el evento, sino que deberá anticiparse en la medida de lo posible, atendiendo a un modelo de predicción de la velocidad de viento. Actualmente y para tal fin, se dispone de datos de predicción del viento “por horas” para municipios, calculados por AEMET a partir de los modelos numéricos de predicción europeos del ECMWF (European Centre for Medium-Range Weather Forecasts). Estas predicciones se ofrecen con un horizonte de hasta 48 horas después de la hora de actualización de dichos modelos (00, 06, 12 y 18 UTC). Para este estudio, la Diputación General de Aragón (DGA) ha facilitado una recopilación de datos horarios de predicción de AEMET a 24 horas, de velocidad y dirección de viento, para el año 2012 y para el municipio de Sabiñánigo. No ha sido posible disponer de predicciones para otras variables meteorológicas, ni para otros periodos de tiempo. En la zona del vertedero de Bailín se ubican una serie de instalaciones pertenecientes a la DGA, destinadas principalmente al estudio y control de la contaminación. Entre ellas, y en la cota menor del vertedero, se encuentra una estación meteorológica fija, gestionada por la empresa GEONICA S.A (ver Figura 6). La estación registra datos cada cinco segundos y calcula distintas medidas de velocidad de viento y dirección asociada, que proporciona en intervalos de 10 minutos. La cantidad de registros es muy amplia y por tanto, los datos en estudio disponibles son muy representativos del régimen de vientos real para la zona. La DGA también ha facilitado los datos recogidos en dicha estación durante el año 2012. 11 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín Como resultado de los estudios anteriores, se comprobará también la adecuación o no de las predicciones de AEMET al comportamiento real del viento. 2.3 ANÁLISIS DE CORRELACIÓN Para finalizar el análisis exploratorio del conjunto de datos, se llevará a cabo un estudio de la correlación entre las covariables de interés, cuyo objetivo es el análisis de las relaciones existentes entre ellas. En el Anexo II se ha incluido una introducción acerca de los conceptos que se utilizan a continuación. Para los pares de variables de velocidad (real y predicha) y dirección (real y predicha), se estudiará el coeficiente de correlación de Pearson y la correlación cruzada; la autocorrelación se estudiará para cada variable por separado. Los análisis se llevarán a cabo mediante gráficos de correlación (correlogramas), herramientas útiles para estudiar la existencia de una correlación retardada (es decir, una correlación mayor con otro instante de tiempo), una autocorrelación, o para analizar la amplitud de los ciclos existentes. 12 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín 3. METODOLOGÍA. CONSTRUCCIÓN DE MODELOS DE REGRESIÓN El fin último de este TFG es realizar un modelo estadístico de predicción de la velocidad del viento que sea capaz de adaptar las predicciones de AEMET al área y periodo de interés. Esta adaptación se ha conseguido mediante la realización de modelos de regresión (MR); en el Anexo II se ha incluido una introducción acerca de los mismos. En este apartado se va a presentar la metodología desarrollada a lo largo del TFG para la construcción de dichos MR. El fin último será encontrar el modelo óptimo de entre todos los construidos, es decir, aquel que ajustándose mejor a la BD en estudio, ofrezca un mejor pronóstico. A lo largo de todo el proceso, y para la implementación de la metodología desarrollada, se utilizará el programa R. Previo al comienzo de la búsqueda del modelo candidato óptimo, se desarrollarán modelos de regresión lineal simple (MRLS), que únicamente incluirán como variable regresora a las predicciones de la velocidad del viento de AEMET. Su objetivo será estudiar la adecuación del modelo de AEMET al régimen real de vientos, de una forma más precisa a la estudiada mediante el análisis exploratorio de datos, y para distintos periodos. Para obtener el MR óptimo se seguirán los pasos que se indican a continuación: 1. Procedimiento de construcción de modelos candidatos, que incluye: a. Estudio y construcción de variables potencialmente útiles. b. Estudio de la transformación de variables. c. Estrategia de construcción de modelos. d. Estudio de la linealidad de ciertas covariables. e. Criterios de bondad de ajuste (en términos de R 2 ajustado). f. Selección de términos (semiautomática y manual) candidatos a formar parte del modelo óptimo. 2. Selección de los “mejores” modelos candidatos una vez desarrollados, atendiendo a distintos criterios. 3. Estudio de la adecuación de los modelos seleccionados mediante el análisis de residuos. 4. Estudio del funcionamiento de los modelos finalmente seleccionados para la realización de pronósticos: a. Estudio gráfico del pronóstico, para distintos umbrales de velocidad. b. Estudio del funcionamiento mediante tablas de contingencia para datos reales y de predicción instantáneos, y estudio de tablas con datos de acumulación horaria para distintos umbrales de velocidad de viento. 13 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín 5. Selección final del modelo óptimo. Desarrollo de un procedimiento en R capaz de obtener pronósticos del modelo óptimo seleccionado a nuevos datos de predicción de AEMET. 3.1 PROCEDIMIENTO DE CONSTRUCCIÓN DE MODELOS CANDIDATOS 3.1.1 ESTUDIO DE VARIABLES POTENCIALMENTE ÚTILES. CONSTRUCCIÓN DE LA BASE DE DATOS En este apartado, se va a explicar cómo desarrollar la BD “final” que contendrá las covariables usadas en la construcción de los MR. A partir de ellas, también podrán considerarse interacciones. Se ha llamado términos al conjunto de covariables e interacciones que formarán parte de distintos modelos candidatos. Se llamará variable o covariable a las variables “aisladas” (es decir, que no forme parte de una interacción). Las covariables que compondrán la BD final son las siguientes:  Covariables que expresen la predicción de la velocidad del viento.  Covariables que expresen la predicción de la dirección del viento.  Covariables que expresen el ciclo diario de brisas.  Covariables que expresen la estacionalidad de la serie de datos. Se tendrán en cuenta las características que las anteriores covariables deban poseer para que puedan ser introducidas al modelo, siendo capaces de reproducir de forma adecuada el comportamiento del fenómeno que con ellas se pretende estudiar. Para la construcción de la BD final, se introducirán en primer lugar los datos de predicción de velocidad de viento de AEMET (instantáneos), así como datos construidos en una ventana temporal alrededor del instante de interés estudiado (t 0 -4, t 0 +4). A dichos datos se les denotará como: -  a la predicción instantánea de velocidad de viento de AEMET, para t 0 . - 1,2,3,4 a los datos de predicción de la velocidad del viento con un desfase horario entre t 0 -1 y t 0 -4 horas (“retardos”, ver Figura 7 izda.). - .1,.2,.3,.4 a la predicción con un desfase horario entre t 0 +1 y t 0 +4 horas (“adelantos”, ver Figura 7 dcha.). Se procederá del mismo modo para los datos de predicción de la dirección, teniendo en cuenta que es una variable circular y viene dada en grados. Por tanto, se trabajará con ella en forma del seno y coseno del primer armónico, aplicando las siguientes expresiones, donde  expresa la dirección predicha en el instante t 0 .  =   2 ;  =   2 Para los armónicos de dirección, se utilizará el mismo código de retardos y adelantos mostrado para las covariables de velocidad. 14 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín Figura 7. Muestra de la construcción de las covariables de predicción de la velocidad del viento alrededor de la ventana temporal (t 0 -4, t 0 +4). A la izda. se muestra el desfase horario existente entre ,1,2,3,4; a la dcha. el desfase entre .1,.2,.3,.4 . En ambas figuras, NA indica que no existen datos. También se utilizarán armónicos para expresar la estacionalidad y los ciclos diarios existentes. Para el ciclo diario se calcularán las siguientes fórmulas, donde ℎ es una componente de la BD que contiene los datos de la hora en punto para cada día, correspondientes a cada observación. .ℎ=   2 ; .ℎ=   2 Por su parte, los términos que representan la estacionalidad se calcularán mediante la siguiente fórmula, donde .ñ nos indica el día del año en que se ha registrado cada dato. Esta covariable tendrá en cuenta desde el día 1 de mayo (que corresponde al día 123 del año 2012) hasta el día 30 de septiembre. . = #.ñ  2 ; . = #.ñ  2 Con todas las covariables anteriores (31 en total, de 4 tipos diferentes) se construirá la BD final que se usará en la construcción de los MR (ver Tabla 1). Tipo de covariable Nº cov. Código en BD Predicción de velocidad de viento 9 instantánea 1  retardos 4  1 ,  2 ,  3 ,  4 adelantos 4  .  1 ,  .  2 ,  .  3 ,  .  4 Predicción de dirección de viento 18 instantánea 2  ,  retardos 8  1 ,  1 ,  2 ,  2 ,  3 ,  3 ,  4 ,  4 adelantos 8  .  1 ,  .  1 ,  .  2 ,  .  2 ,  .  3 ,  .  3 ,  .  4 ,  .  4 Ciclo diario 2  . ℎ  ,  . ℎ  Ciclo anual 2  .  ,  .  TOTAL 31 Tabla 1. Base de datos final, usada para el desarrollo de MR candidatos a ser el modelo de predicción óptimo. 15 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín 3.1.2 ESTUDIO DE UNA POSIBLE TRANSFORMACIÓN DE VARIABLES El estudio de la transformación de la respuesta será necesario cuando los residuos de los MR en estudio no cumplan las hipótesis de distribución normal y varianza constante (ver apartado 3.3). Para ello, se estudiará la transformación Box Cox, procedimiento estadístico que se resume en el Anexo II. Se considerará también una transformación de tipo polinómico para las covariables regresoras de predicción de la velocidad de viento, las únicas susceptibles de precisar de dicha transformación para explicar su relación con la variable respuesta. 3.1.3 ESTRATEGIA DE CONSTRUCCIÓN DE MODELOS El elevado número de potenciales covariables hace inviable la construcción exhaustiva de todos los modelos posibles. Por ello, se ha diseñado una estrategia que ha consistido en el desarrollo de modelos sucesivos, desde los más sencillos posibles (con tan solo un tipo de covariables y sin interacciones) hasta los más complejos, con todos los tipos de covariables, transformaciones polinómicas e interacciones de todo tipo.  Desarrollo de modelos de regresión lineal sin interacciones En primer lugar se desarrollarán modelos sucesivos, a los cuales se irán añadiendo, en cada paso, covariables de diferente tipo. En el proceso se estudiará cuáles de dichas covariables son más significativas para explicar la respuesta. Las covariables añadidas serán, en este orden, covariables de predicción de la velocidad del viento, armónicos de predicción de la dirección del viento, armónicos de estacionalidad y armónicos de ciclo diario. A su vez, las covariables se añadirán mediante dos métodos diferentes. El método 1 consistirá en añadir un nuevo tipo de covariable al “mejor” modelo hallado en el paso anterior; el método 2 estudiará directamente el modelo que incluya todas las covariables potenciales, en un único paso. Para hallar el “mejor” modelo en ambos casos, se seguirán los criterios de bondad de ajuste y de selección de variables, posteriormente explicados en el apartados 3.1.4. Mediante cada método se obtendrá un “mejor modelo”, susceptible de contener los cuatro tipos de covariables. Se realizará entonces un estudio de la linealidad de los “mejores modelos”, es decir, de la necesidad de expresiones polinómicas para las covariables de predicción de la velocidad del viento. El proceso consistirá en realizar distintos modelos anidados (es decir, donde los predictores de los modelos más sencillos sean un subconjunto de los que forman los modelos más complejos) con la función gam, generalized additive models (Hastie and Tibshirani, 1991), y comparar dichos modelos con un test ANOVA (Fox and Weisberg, 2011), comprobando cuál de ellos resulta más significativo. En el Anexo II se ha recogido una breve introducción acerca del test ANOVA. Tras el estudio de la linealidad, se desarrollarán nuevos modelos en los cuales se incluirán los términos polinómicos que hayan resultado significativos. 16 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín  Desarrollo de modelos de regresión con interacciones Una vez que se hayan obtenido los “mejores” modelos posibles sin interacciones, lineales y polinómicos, y con el fin de aumentar la bondad de ajuste del modelo, se estudiará la introducción de interacciones a los modelos candidatos, para encontrar aquellas que puedan resultar significativas. La introducción de interacciones a los modelos podrá realizarse mediante bucles, lo cual permitirá estudiar un espectro de combinaciones de covariables e interacciones mucho mayor que mediante la introducción manual. Las interacciones que se ensayarán serán las que a priori, y dados los análisis exploratorios anteriores, se crea que resultarán más significativas. La búsqueda mediante bucles será de dos tipos: - Bucle 1  Considerará covariables e interacciones guardadas en un vector, una en cada repetición del bucle, y desarrollará un modelo con cada una mediante un procedimiento “paso a paso” (stepwise, ver apartado 3.1.4), partiendo de uno de los “mejores” modelos obtenidos en el apartado anterior. En cada repetición, seleccionará y guardará en otro vector todos aquellos términos que hayan resultado significativos. Una vez el bucle termine, se ajustará manualmente un MR que incluya todos los términos guardados. - Bucle 2  Considerará covariables e interacciones guardadas en el mismo vector anterior, y desarrollará un modelo con ellas. La diferencia es que en este caso, el modelo se “actualizará”, es decir, en cada repetición del bucle el modelo añadirá nuevos términos y eliminará aquellos que no resulten significativos, también en un procedimiento de tipo stepwise. Por tanto, el resultado del bucle será directamente un MR candidato, resultado de aplicar todos los términos en estudio, en todas las repeticiones. A lo largo del desarrollo de los modelos es importante tener en cuenta que habrá que seleccionar, tanto términos a incluir en el modelo, como modelos propiamente dichos, mediante criterios de selección de variables y de bondad de ajuste respectivamente, los cuales se estudian a continuación. 3.1.4 CRITERIOS DE BONDAD DE AJUSTE Y SELECCIÓN DE VARIABLES Los criterios de bondad de ajuste permiten medir la capacidad del modelo para explicar el comportamiento de la variable respuesta. De entre los existentes, se ha seleccionado el R 2 ajustado. En el Anexo II se puede encontrar más información sobre el mismo. Anteriormente, se ha explicado cómo construir la BD final, constituida por 31 covariables, que junto con las interacciones construidas entre ellas, serán candidatas a formar parte del MR óptimo que se busca. Sin embargo, probablemente solo un subconjunto de ellas estará relacionado con la respuesta. Por tanto, se debe realizar un proceso de selección, cuyo fin será encontrar la “mejor ecuación de regresión” (Montgomery et al., 2012), es decir, un MR con suficientes términos como para capturar el comportamiento de la variable respuesta, pero que 17 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín sea fácil de usar e interpretar y por tanto, que tenga el menor número posible de ellos. La selección se puede llevar a cabo de forma exhaustiva o automática. La aproximación exhaustiva es la más directa a la hora de decidir qué términos son importantes para el modelo, y consiste en estudiar todas las combinaciones de modelos posibles para encontrar aquellas más predictivas. Este estudio solo resulta abordable si el número de términos en estudio es pequeño, ya que las posibles combinaciones con  parámetros es del orden de 2 $ (James et al., 2013). Por tanto, cuando se dispone de un número elevado de términos, se necesita de una aproximación semiautomática. En el Anexo II se puede encontrar una introducción acerca del método semiautomático stepwise de selección mixta, utilizado en este trabajo. Cuando se lleva a cabo la selección mixta puede ocurrir que el modelo resultante contenga covariables no significativas, porque formen parte de interacciones que sí lo sean. Esto ocurre por el llamado principio de jerarquía, que expone que si se incluye una interacción en el modelo, también deben incluirse las covariables que la forman, para que dicho modelo se pueda interpretar (James et al., 2013). Sin embargo, se ha comprobado que los modelos candidatos son demasiado complejos en términos de interpretabilidad. Además, se ha establecido como criterio que todos los parámetros regresores incluidos en el modelo deberán ser significativos para el mismo, es decir, tener un p-valor asociado a sus parámetros inferior a 0.05. Esto supone que tras la aplicación de la selección mixta, se deberá realizar una limpieza manual del modelo de términos no significativos para el mismo. Esta actualización del modelo se debe hacer “paso por paso”, es decir, se debe eliminar (o añadir) un único término cada vez, estudiando cómo afecta este cambio a la significación, tanto del modelo como de los términos que lo forman, hasta obtener un modelo en el cual todos los términos sean significativos. En caso de que existan varios modelos anidados, con un ajuste similar y con distinto número de variables, se estudiará su adecuación mediante su comparación con un test ANOVA. 3.2 SELECCIÓN DE MODELOS CANDIDATOS Una vez se haya desarrollado un conjunto de modelos candidatos, sin y con interacciones, se hará una selección de los “mejores” atendiendo a los siguientes criterios:  La bondad de ajuste de cada modelo, expresada por su R 2 ajustado y por su desviación típica residual.  La composición de los mismos, es decir, el número de términos que lo forman, el grado de significación de dichos términos a la hora de predecir la respuesta, la existencia de parámetros no deseables o el número de interacciones triples.  La desviación típica residual de validación cruzada, técnica posteriormente introducida.  La bondad de ajuste de cada modelo a cada mes por separado y en cada una de las 24 horas. 18 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín El proceso de selección de los modelos candidatos consiste en hallar un equilibrio entre los criterios anteriores; los criterios que primarán serán la bondad de ajuste y el resultado de la validación cruzada, en términos de desviación típica residual. En especial, si ambos criterios tienen similar resultado para los modelos en estudio, será importante atender a su composición, teniendo preferencia los modelos más sencillos, evitando con ello posibles sobreajustes. El estudio de la composición también aportará información sobre posibles deficiencias en el comportamiento del modelo, que deberán subsanarse. Se atenderá asimismo al grado de significación de los términos, prefiriendo modelos cuyos términos sean muy significativos a la hora de explicar la respuesta (es decir, cuyo p-valor sea inferior a 0.01). Por último, la bondad de ajuste de cada modelo a cada mes por separado, indicará si el modelo funciona correctamente para todo el periodo, o si por el contrario será necesario el desarrollo de un modelo adaptado a algún mes o periodo del día en particular. 3.3 ESTUDIO DE LA ADECUACIÓN DE LOS MODELOS SELECCIONADOS La adecuación de los modelos seleccionados se estudiará mediante el análisis de sus residuos, que permiten comprobar si el modelo tiene realmente capacidad predictiva, es decir, si es capaz de recoger el comportamiento de la variable respuesta en estudio. Para considerar que la estimación de los parámetros de los modelos es adecuada, se deben suponer, al menos, las siguientes hipótesis: - Los errores tienen media cero. - Los errores tienen varianza constante (situación de homocedasticidad). - Los errores se distribuyen según una normal. El cumplimiento de las hipótesis anteriores es condición básica para considerar como apropiado el test de hipótesis, que determina si existe relación lineal entre la respuesta y un subconjunto de variables regresoras, siendo: %  :' ( ='  =⋯=' * =0 % ( :' , ≠0 , para al menos un valor de j Así, el rechazo de la hipótesis nula implica que al menos una de las variables regresoras contribuye al modelo de forma significativa. El análisis de los residuos se realizará mediante gráficos de diagnóstico que permitirán el estudio de su linealidad, homocedasticidad, distribución y puntos influyentes. La distribución se estudiará también mediante histogramas. Todos los gráficos anteriores vienen explicados en el Anexo II. En caso de encontrarse deficiencias en la normalidad de los residuos de los modelos candidatos seleccionados, se llevarán a cabo los siguientes estudios: 19 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín - Transformación Box Cox de la respuesta. - Gráfico secuencial de residuos. - Gráfico de residuos sobre variables regresoras incluidas en el modelo (velocidad y dirección del viento). - Gráficos de residuos en forma de diagramas de caja, para meses y horas. - Gráficos residuales parciales, obtenidos en R mediante la función crPlots (Fox and Weisberg, 2011), que evalúan el efecto de cada predictor sobre la variable respuesta. Por tanto, permiten estudiar deficiencias en el modelo (p.e., necesidad de transformaciones polinómicas). - Gráficos de variable añadida, obtenidos en R mediante la función avPlots (Fox and Weisberg, 2011), en los cuales se estudiará la pendiente hallada en cada gráfico, pues ésta representa el coeficiente parcial de regresión para el parámetro en estudio. 3.4 ESTUDIO DE LA CAPACIDAD PREDICTIVA DEL MODELO Anteriormente, se han citado los criterios de selección de modelos estudiados. Uno de ellos era la desviación típica residual de validación cruzada. Se llama validación a la comprobación de la capacidad predictiva del modelo, es decir, si los modelos ajustados para la BD en estudio, funcionarían también para un nuevo conjunto de observaciones, o si por el contrario el modelo está sobreajustado y por tanto, solo es válido para la BD que se ha usado en su ajuste. Antes de estudiar cómo realizar dicha validación, se deben tener en cuenta dos conceptos importantes: el error de validación (test error) y el error de entrenamiento (training error). Ambos se han explicado en el Anexo II. Cuando se quiere validar un MR existen dos alternativas. La primera consiste en estudiar cómo se comporta el modelo en una BD que no se haya utilizado en su ajuste, con lo que se obtiene directamente el error de validación. La alternativa, cuando no existe una nueva BD, es aplicar herramientas estadísticas de remuestreo (resampling) que permiten validar un modelo usando las observaciones previamente utilizadas para su ajuste. De los métodos de remuestreo existentes, se ha seleccionado el método de validación cruzada para k subconjuntos, que se da en dos pasos: 1. División de la BD en k subconjuntos, de forma consecutiva (ver ejemplo en Figura 8), o aleatoria (ver ejemplo en Figura 9). 2. Reajuste del modelo, en k iteraciones; en cada iteración se toma un subconjunto de prueba, y un subconjunto de entrenamiento. En el Anexo II se ha incluido una explicación más desarrollada acerca del proceso de validación cruzada para k subconjuntos. La validación cruzada se llevará a cabo en R mediante el paquete cvTools (Alfons, 2012). 20 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín Figura 8. Ejemplo gráfico del método consecutivo de validación cruzada, en 5 iteraciones (It.). Figura 9. Ejemplo gráfico del método aleatorio de validación cruzada, en 5 iteraciones (It.). 3.5 ESTUDIO DEL FUNCIONAMIENTO OPERATIVO DE LOS MODELOS SELECCIONADOS. REALIZACIÓN DE PRONÓSTICOS Una vez se ha estudiado la adecuación de los modelos seleccionados, comprobando que se cumplen las hipótesis sobre los residuos y que los modelos se ajustan a “nuevas” BD, y previo a la selección final del modelo óptimo, se estudiará el funcionamiento real de los modelos. Resulta muy complicado que un modelo de tipo meteorológico acierte en sus predicciones de forma instantánea. Es decir, que se cumplan los valores que el modelo ha predicho para el momento exacto para el cual los ha predicho. Por tanto, no interesa tanto estudiar las predicciones instantáneas de los modelos desarrollados, sino su comportamiento en términos de probabilidad de ocurrencia del fenómeno en estudio (pronósticos). En este caso, se estudiará la probabilidad de que el valor de velocidad de viento 27 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín En las rosas de vientos, la escala de velocidad de viento estudiada ha ido desde 0 y hasta el percentil 96.1 correspondiente a cada BD, en intervalos de 5 km/h para los datos reales y de 3.24 km/h para los predichos; de esta forma ambas tienen el mismo número de intervalos y el mismo código de color y son comparables gráficamente. La comparación por intervalos se ha hecho para las llamadas situaciones 1 y 2. Recordemos que la situación 1 estudiaba los datos reales para cada intervalo y sus predicciones instantáneas correspondientes, mientras la situación 2 estudiaba los datos predichos pertenecientes a cada intervalo y las observaciones instantáneas correspondientes. En la Figura 13 se recogen las rosas de vientos elaboradas para el estudio de la situación 1, y para los intervalos de velocidad del viento correspondientes a intensidades altas y medias. Para el intervalo de intensidades elevadas de viento, los datos reales muestran una clara componente WNW (ver Figura 13a), al igual que sus predicciones instantáneas, las cuales en general muestran velocidades elevadas (ver Figura 13b). Figura 13a Figura 13b Figura 13c Figura 13d Figura 13. Rosas de vientos correspondientes al estudio de la situación 1 para intervalos de intensidades elevadas y medias de viento. La Figura 13a. corresponde a datos observados en el intervalo de intensidad de viento elevada y la Figura 13b muestra sus datos de predicción correspondientes. La Figura 13c muestra datos observados en el intervalo de intensidad de viento media, mientras la Figura 13d. recoge sus datos de predicción correspondientes. Para intensidades de viento medias, los datos reales siguen mostrando una componente WNW predominante, aunque aparecen nuevas componentes, en especial W y ESE (ver Figura 13c). Por 28 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín su parte, los datos de predicción correspondientes muestran una componente W predominante, además de muchas otras que no se observan en el régimen real de vientos (ver Figura 13d). Se predicen velocidades más altas de las que realmente se observan. Para situaciones de calma, el régimen real de vientos cambia radicalmente de componente principal, pasando a ser SSE (ver Figura A. 3a). Sin embargo, se observan velocidades muy elevadas para la predicción instantánea correspondiente (ver Figura A. 3b). La componente principal de la predicción de la dirección también es diferente en este caso, destacando la componente E. En el estudio de la situación 2, recogido en la Figura 14, se observa que para intensidades elevadas de viento, las velocidades predichas muestran una componente E mayoritaria, siendo importante también la componente WNW (ver Figura 14a). Los datos reales correspondientes muestran, al igual que en la situación 1 para el mismo intervalo, una marcada componente WNW, aunque con un mayor rango de velocidades (Figura 14b). No se observa la componente E en las velocidades reales. Figura 14a Figura 14b. Figura 14c. Figura 14d. Figura 14. Como Figura 13, para estudio de rosas de vientos correspondientes a situación 2. En este caso Figuras 14 a y c muestran predicciones, y Figuras 14 b y d observaciones correspondientes Para intensidades medias de viento, la dirección preferente del modelo de predicciones de AEMET es E de nuevo (Figura 14c), mientras que los datos reales correspondientes muestran una clara dirección WNW y SSE (Figura 14d). Se observa que el modelo de predicciones de AEMET, aparte de no poder capturar el comportamiento del régimen real de vientos en términos 29 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín de dirección, tampoco predice las velocidades elevadas que se observan en realidad para los mismos instantes de tiempo. Para el intervalo de calmas, se observan predicciones para prácticamente todas las direcciones contenidas en la rosa de vientos (Figura A. 4a), mientras las observaciones correspondientes instantáneas muestran una dirección preferente SSE (Figura A. 4b), al igual que ocurría en la situación 1 para el mismo intervalo de vientos. También se observa que cuando el modelo de AEMET predice situaciones de calma, las observaciones reales de velocidad de viento no se encuentran necesariamente en dicho intervalo, sino que muestran un rango mucho más amplio, llegando a observarse velocidades de viento muy elevadas que como decimos, no quedan reflejadas en las predicciones de AEMET. En general, se han encontrado grandes diferencias entre las rosas de vientos de datos observados y datos instantáneos predichos por AEMET, para cada intervalo y situación en estudio. Por tanto, de nuevo se puede concluir que el modelo de predicciones de AEMET no se adecua correctamente al régimen real de vientos. ESTUDIO DE LA ESTACIONALIDAD Y DEL CICLO DIARIO DEL RÉGIMEN DE VIENTOS Siguiendo con la metodología desarrollada en el apartado 2.2, se ha realizado un estudio de la estacionalidad y del ciclo diario de las BD en estudio. Se ha descartado profundizar en el estudio de la tendencia de la BD.max2, ya que un análisis previo de tipo descomposición de serie de datos, no ha arrojado resultados significativos sobre la existencia de dicha tendencia (ver Figura A. 43 en Anexo III, sobre estudios adicionales).  Estudio del ciclo diario El comportamiento del ciclo diario se ha estudiado mediante un gráfico de variación horaria de las velocidades medias durante el periodo en estudio, reales y predichas, y mediante la comparación de las rosas de los vientos para los horarios diurnos (entre 8 y 19 horas) y nocturnos (entre 20 y 7 horas). La Figura 15 muestra el gráfico de variación horaria de las velocidades medias. En él se observan claras diferencias entre datos reales y predichos: el régimen real de velocidad muestra la existencia de un ciclo diario, el cual tiene un máximo alrededor de las 15 horas (máximo que también se observa para el modelo de predicciones) y disminuye, siguiendo una forma sinusoidal hasta alcanzar un mínimo a las 8 de la mañana. Se observa un ciclo menos marcado para el modelo de predicciones de AEMET, manteniéndose el nivel medio más homogéneo, con un rango en la velocidad media (entre los valores mínimo y máximo) de hasta 6 km/h; este rango aumenta hasta los 16 km/h para los datos reales (ver Figura 15). 30 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín Figura 15. Estudio de la variación de las velocidades medias real y predicha a lo largo del día. Las líneas discontinuas indican los puntos de cambio del ciclo diario existente en el régimen real de velocidades. Para el horario diurno, las rosas de los vientos muestran cómo el régimen real de vientos tiene una clara componente WNW y una ligera componente ESE (ver Figura 16 izda.), con velocidades más o menos elevadas. Por su parte, el modelo de predicciones de AEMET muestra un régimen de vientos con predominancia de la componente W (ver Figura 16 dcha.), así como una distribución similar de frecuencias de velocidad de viento. Figura 16. Rosas de vientos realizadas con datos correspondientes al horario diurno. A la izda. Se muestra la rosa para BD.max2 y a la dcha., la correspondiente a BD.p, cada una con una escala adecuada. La distribución de frecuencias de BD.max2 para el horario diurno es de tipo campana de Gauss (ver Figura A. 5a), con una frecuencia máxima del 28% para el intervalo de 15-20 km/h. Las velocidades superiores a 20 km/h suponen un 34.5% en total, siendo superiores a 30 km/h un 6.8%. La distribución de las frecuencias BD.p para el mismo horario (ver Figura A. 5.b) tiene una forma similar, pero se encuentra desplazada hacia la izquierda, mostrando el suavizado del 0 5 10 15 20 0 5 10 15 20 25 30 Variación diaria v. media Hora v.media (km/h) v. media real v. media predicha 31 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín predicciones de AEMET con respecto a los datos observados. En este caso, los porcentajes máximos de frecuencias están concentrados en los intervalos entre 5-10 km/h (34.2%) y 10-15 km/h (35.7%). El intervalo de velocidades de 15-20 km/h agrupan un 11.4% de los datos, mientras que solo un 3.7% de los datos son superiores a 20 km/h, y un 0.7% a 30 km/h. Es decir, para velocidades superiores a 20 km/h, las frecuencias del modelo de predicciones son aproximadamente una décima parte de las realmente observadas (ver Tabla 2). HORARIO DIURNO FRECUENCIA HORARIO NOCTURNO FRECUENCIA Intervalo velocidad (km/h) Observado Predicho Intervalo velocidad (km/h) Observado Predicho 5 a 10 - 34.2% 5 a 10 51.7% 27.1% 10 a 15 - 35.7% 10 a 15 24.7% 25.2% 15 a 20 28% 11.4% - - - > 20 34.5% 3.7% >15 14.7% 14.4% > 30 6.8% 0.7% > 30 1% - Tabla 2. Estudio del porcentaje de frecuencias de acumulación de datos de velocidad del viento para cada intervalo de viento estudiado, horarios diurno (izda.) y nocturno (dcha.). El símbolo “-” indica que los datos para ese intervalo no son relevantes para el estudio. El comportamiento real del viento durante la noche cambia radicalmente, siendo la dirección preferente SSE, y registrando bajas velocidades en comparación con las observadas durante el día (ver Figura 17 izda.). El modelo de predicciones de AEMET muestra un régimen de vientos con una marcada componente E (Figura 17 dcha.), y rangos de velocidades mucho mayores de los que se dan realmente. Figura 17. Rosas de vientos realizadas con datos correspondientes al horario nocturno. A la izda. Se recoge la rosa de vientos para BD.max2; a la dcha., la correspondiente a BD.p. Para el horario nocturno, la distribución de frecuencias de intervalos de velocidad de viento cambia respecto al diurno, siendo en este caso distinta para ambas BD (ver Figura A. 5 c y d). Para BD.max2 (ver Tabla 2), un 51.7% de las observaciones pertenecen al intervalo de 5-10 km/h, un 24.7% de 10-15 km/h, y un 14.7% entre 15-30 km/h; únicamente un 1% de las observaciones superan los 30 km/h en comparación con el 6.8% obtenido para el horario diurno. Las predicciones instantáneas muestran una distribución más homogénea de frecuencias (ver Figura A. 5.d), 27.1% para 5-10 km/h, 25.2% para 10-15 km/h y un total de 14.4% para velocidades superiores a 15 km/h (ver Tabla 2). 32 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín En resumen, para el horario nocturno y en cuanto a la distribución de la velocidad, el modelo de predicciones de AEMET se comporta correctamente para velocidades superiores a 10 km/h, probablemente dada la menor variabilidad de las observaciones registradas durante la noche; para el horario diurno y para velocidades observadas superiores a 20 km/h, el modelo de predicciones de AEMET recoge únicamente un ≈10% de ellas.  Estudio de la estacionalidad La estacionalidad se ha estudiado mediante gráficos de variación de la velocidad media a lo largo del periodo en estudio, y de dicha variación para distintos cuartiles. En el gráfico de variación de la velocidad media mensual, recogido en la Figura 18, se observa que la evolución que siguen las medias de los datos de velocidad real y predicha es muy diferente: la velocidad media real varía de forma estable a lo largo del periodo, siguiendo una curva que tiene un máximo en agosto para posteriormente descender, mientras que la velocidad media predicha para cada mes distingue dos mínimos relativos en julio y septiembre. Figura 18. Estudio de la variación de los datos medios de velocidad real y predicha para cada mes del periodo entre mayo y octubre. Esta diferencia aparente entre la velocidad de viento real y predicha, también aparece en el estudio de los distintos cuartiles para ambas BD. En la Figura 19 izda. se observa que la tendencia que siguen los percentiles 0.25 y 0.50 de las velocidades medias es bastante similar, aunque no lo es para el percentil 0.75. Por el contrario, si estudiamos la tendencia de los cuartiles para la velocidad predicha (ver Figura 19 dcha.) todos ellos muestran una evolución casi paralela, simetría que no se observa en el gráfico de cuartiles de velocidad real. 8 10 12 14 16 Velocidad media real y predicha (km/h) Meses Velocidad media (km/h) mayo junio julio agosto septiembre octubre v. media real v. media predicha 33 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín Figura 19. Estudio de las velocidades real (izda.) y predicha (dcha) correspondientes a los 3 cuartiles (0.25, 0.5 y 0.75). 4.1.3 ANÁLISIS DE CORRELACIÓN Como último paso del análisis exploratorio de datos, se ha estudiado la correlación entre las variables reales y predichas según lo especificado en el apartado 2.3. El fin de este análisis es valorar la existencia de una relación lineal entre las variables en estudio, así como comprobar si existe algún tipo de retardo en dicha relación, es decir, si los datos de predicción podrían estar desfasados por algún motivo.  Velocidad del viento real y predicha Para la pareja de covariables formada por la velocidad del viento real y la predicha, el coeficiente de correlación de Pearson hallado es de 0.41 y es significativa a nivel 0.05. Atendiendo a la Figura 20, que muestra el gráfico de correlación cruzada entre ambas covariables, vemos que el valor máximo de correlación corresponde al instante 0 (lag=0 en el gráfico), estableciendo que no hay un desfase sistemático entre las covariables a lo largo del tiempo. El comportamiento de la correlación entre ambas covariables sigue una forma sinusoidal, con un ciclo completo de aproximadamente 24 horas, tras el cual la correlación cruzada es de nuevo máxima y positiva (ver línea verde en el gráfico), dividido en ciclos de 12 horas en los cuales la correlación pasa a tener signo negativo (línea discontinua roja en el gráfico). Esta situación expresa el ciclo diario de brisas existente (la diferencia en la velocidad del viento entre el día y la noche, ya observada en análisis anteriores), que se refleja en ambas series de datos, aunque menos nítidamente en la predicha (ver Figura 21). Los gráficos de autocorrelación de las variables muestran cómo la velocidad real del viento cambia el signo de su autocorrelación cada 12 horas, retomando el signo en ciclos de 24 horas (ver Figura 21.izda). Sin embargo, las predicciones de AEMET de velocidad del viento muestran un ciclo completo cada 24 horas, pero no un cambio en el signo de su autocorrelación (ver Figura 21.dcha). 0 5 10 15 20 25 Velocidad real (km/h) 0 5 10 15 20 25 Velocidad predicha (km/h) 34 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín Figura 20. Gráfico de correlación cruzada entre velocidad del viento real y predicha. Figura 21. Conjunto de correlogramas para variables de velocidad del viento real (v.viento, izda.) y predicha (p.v.viento, dcha.).  Dirección del viento real y predicha Se ha realizado el mismo estudio comparando las variables de dirección del viento real y su predicción correspondiente. En este caso, aunque se trata de una variable circular, la rosa de vientos reales ha indicado la ausencia de viento de componente Norte (ver Figura 13), por lo que se ha concluido que se puede trabajar con la correlación como una herramienta descriptiva de la relación. El coeficiente de correlación de Pearson hallado para ambas variables es de tan solo un 0.26. El gráfico de correlación cruzada de ambas variables (ver Figura 22), muestra un máximo en lag=-2 (≈0.35), es decir, con un desfase horario de 2 horas. Al igual que ocurría con la velocidad -30 -20 -10 0 10 20 30 -0.1 0.0 0.1 0.2 0.3 0.4 Lag ACF Correlación entre v.viento y p.v.viento 35 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín del viento, el gráfico tiene una clara forma sinusoidal que cambia la tendencia de la correlación en ciclos de aproximadamente 12 horas, retomando el signo positivo en ciclos de 24 horas. Figura 22. Gráfico de correlación cruzada entre dirección del viento real y predicción de la dirección del viento . Si se comparan los autocorrelogramas para ambas covariables, se observa que la dirección real muestra un patrón de giro del signo de la autocorrelación cada 12 horas (ver Figura 23.izda), mientras que las predicciones son capaces de mostrar autocorrelación nula en el retardo 12, pero no negativa (ver Figura 23 dcha.). Figura 23. Conjunto de gráficos de autocorrelación para variables de dirección real del viento (dir.viento, izda.) y predicha (p.dir.viento, dcha.) Los resultados obtenidos en el análisis de correlación de nuevo muestran la existencia de un régimen de brisas diferenciado para el horario diurno y nocturno, el cual el modelo de predicciones no es capaz de reflejar satisfactoriamente. -30 -20 -10 0 10 20 30 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 Lag ACF Correlación entre dir.viento y p.dir.viento 36 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín 4.1.4 PRINCIPALES CONCLUSIONES DEL ANÁLISIS EXPLORATORIO Los análisis exploratorios realizados hasta el momento sugieren descartar el uso directo del modelo de predicciones de AEMET, por las carencias que dicho modelo ha mostrado para el área del vertedero de Bailín y el periodo de verano. Será necesario por tanto, desarrollar un nuevo modelo estadístico de predicciones de la velocidad del viento. Dicho modelo deberá ser capaz de reescalar las predicciones de AEMET, incluyendo componentes que reflejen la estacionalidad y el ciclo diario existente en el régimen real de vientos, características ausentes en las predicciones de AEMET. Además, el nuevo modelo deberá incluir términos de interacción entre las covariables anteriores, dada la relación observada entre por ejemplo, la velocidad y la dirección del viento, la velocidad del viento y el ciclo diario de brisas, o la velocidad del viento y la estacionalidad. En el análisis de la correlación se ha observado una situación en la cual ésta es mayor para una diferencia temporal de -2 horas (entre las direcciones de viento reales y predichas). Descartando que esta diferencia pudiese deberse a que los datos de predicción se diesen de forma decalada, surgió la idea de crear “nuevas” covariables a incluir en el modelo, a partir de los datos de predicción de AEMET (tanto de velocidad como de dirección), en una ventana temporal de 4 horas alrededor del instante de interés estudiado (t 0 -4,t 0 +4), ya que en dicho intervalo la correlación cruzada se mantiene superior a 0.2, tanto para la velocidad como para la dirección. Por último, en este nuevo modelo estadístico de predicciones deberá estudiarse la conveniencia de la inclusión del mes de octubre, ya que se ha observado que posee un régimen de vientos más homogéneo que el resto de periodo en estudio. 4.2 DESARROLLO DE MODELOS DE REGRESIÓN Según lo expuesto en el apartado 3.1 de la metodología, el tipo de modelos seleccionados para realizar la adaptación (downscaling) de las predicciones de AEMET al periodo y área de interés, han sido los modelos de regresión (MR). Los principales pasos seguidos para la construcción de los modelos han sido los siguientes: 1. Estudio de la capacidad para distintos periodos, del modelo de regresión lineal simple (MRLS) construido a partir de las predicciones de AEMET. 2. Desarrollo de MR candidatos a ser el modelo estadístico óptimo buscado. 3. Selección de modelos candidatos. 4. Adecuación de los modelos seleccionados. 5. Estudio del funcionamiento operativo de los modelos seleccionados. Pronóstico. 43 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín Por tanto, se ha decidido trabajar también con los dos modelos que incluyen polinomios. MODELO Estudio nolinealidad R 2 ajustado Nº predictores modelo.original * 0.13 1 ensayo1 * 0.16 3 ensayo2 * 0.34 11 ensayo3 * 0.33 12 ensayo4 0.47 10 ensayo4.poly 0.48 16 ensayo4.1 0.47 14 ensayo4.1.poly 0.49 16 Tabla 4. Resumen de los modelos obtenidos mediante el desarrollo de MR sin interacciones. En general, se puede concluir que para el estudio de los MR sin interacciones, se da una mejora significativa en el R 2 ajustado, llegando éste a alcanzar en torno al 0.48, cuando se añaden: - Covariables ligadas a la dirección del viento, que aumentan la explicación de la variabilidad de la respuesta en un 18%. - Covariables que expresan el ciclo diario existente en el régimen de brisas, que mejoran la explicación de la respuesta en un 13%. - Los términos polinómicos únicamente aumentan la explicación de la variabilidad de la respuesta en un 1%, aunque según los test ANOVA realizados, resultan significativos. Los modelos candidatos que se utilizarán para seguir desarrollando el modelo “óptimo” son de tipo lineal (ensayo4 y ensayo4.1) y polinómico (ensayo4.poly y ensayo4.1.poly). DESARROLLO DE MODELOS DE REGRESIÓN CON INTERACCIONES Se pretende incrementar el grado de explicación de la respuesta (velocidad observada) considerando la introducción de interacciones entre componentes que intervienen en el MR. El efecto de las interacciones se ha considerado construyendo modelos dirigidos por las covariables existentes y también mediante procedimientos semiautomáticos.  Introducción manual Siguiendo la metodología, y dados los indicios hallados en el análisis exploratorio, las interacciones evaluadas han sido las siguientes: - Interacciones entre covariables de predicción de la velocidad del viento y estacionalidad, ya que el régimen de vientos cambia a lo largo del año. - Interacciones entre covariables de predicción de la velocidad y predicción de la dirección del viento, ya que se ha detectado una relación entre ambas. - Interacciones entre covariables de predicción de la velocidad y ciclo diario, ya que se ha comprobado la existencia de un régimen diario de brisas diferenciado. - Interacciones entre ciclo diario y estacionalidad. 44 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín El test ANOVA muestra como modelos más significativos aquellos a los cuales se introducen todas las interacciones en estudio (ver Figura A.13 para modelos lineales y Figura A.14 para polinómicos). Dichos modelos tienen mayor R 2 ajustado y número de parámetros. Sin embargo, la interacción que produce el mayor aumento en la bondad de ajuste del modelo con el menor aumento en el número de parámetros (la “mejor ecuación de regresión” en este caso), es aquella formada por las covariables de predicción de velocidad de viento y ciclo diario (ver casillas coloreadas en Tabla 5). Este aumento es de un 2.5% con solo 1 parámetro más para los dos modelos lineales, y de un 2.5% para los polinómicos, aunque con un mayor aumento en el número de parámetros (entre 6 y 11, ver Tabla 5). MODELOS LINEALES N.par. R 2 ajustado Interacción añadida MODELOS POLINÓMICOS N.par. R 2 ajustado Modelo original ensayo4 10 0.471 - - Modelo original ensayo4.poly 16 0.485 ensayo4.a 13 0.473 p.v.viento*ciclo.anual ensayo4.poly.a 26 0.492 ensayo4.b 11 0.496 p.v.viento*ciclo.diario ensayo4.poly.b 27 0.511 ensayo4.c 12 0.497 p.v.viento*(ciclo.diario + ciclo.anual) ensayo4.poly.c 32 0.515 ensayo4.d 14 0.499 ensayo4.c + ciclo.diario*ciclo.anual ensayo4.poly.d 34 0.517 Modelo original ensayo4.1 14 0.474 - - Modelo original ensayo4.1.poly 16 0.486 ensayo4.1.a 17 0.477 p.v.viento*ciclo.anual ensayo4.1.poly.a 23 0.493 ensayo4.1.b 15 0.500 p.v.viento*ciclo.diario ensayo4.1.poly.b 22 0.511 ensayo4.1.c 15 0.501 p.v.viento*(ciclo.diario + ciclo.anual) ensayo4.1.poly.c 31 0.517 ensayo4.1.d 19 0.503 ensayo4.1.c + ciclo.diario*ciclo.anual ensayo4.1.poly.d 23 0.512 Tabla 5. Modelos lineales y polinómicos a los que se han introducido interacciones de forma manual. El símbolo (- -) indica que el modelo no tiene ninguna interacción. Se han coloreado los modelos a los que se han introducido interacciones que producen una mejora más significativa.  Introducción mediante procedimientos semiautomáticos El segundo ensayo se basa en la introducción de interacciones a los modelos mediante el uso de los bucles desarrollados en R para tal fin (b1 y b2, ver apartado 3.1.3), que han permitido estudiar una gran cantidad de modelos diferentes de forma automática y rápida. Los modelos obtenidos se han resumido mediante el nombre del modelo de partida seguido del bucle utilizado en cada caso (b1, b2). Las interacciones (hasta de orden 3) se han estudiado entre las covariables asociadas a velocidad y dirección predichas y armónicos de los ciclos diarios y estacionales; esto corresponde a aproximadamente 5000 términos potenciales distintos. Si se comparan los modelos obtenidos mediante b1 y b2, tanto a partir de modelos lineales como polinómicos (ver Tabla 6), se observa cómo el uso b2 obtiene modelos con menor número de parámetros y similar bondad de ajuste. 45 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín La realización de un test ANOVA, muestra resultados más significativos para el ensayo4.1.b2 de entre los modelos lineales (ver Figura A.15) y para el ensayo4.poly.b2 de entre los polinómicos (ver Figura A.16), aunque en este caso, no se ha podido comparar con el ensayo4.1.poly.b2 por no estar anidados. MODELOS LINEALES MODELOS POLINÓMICOS BUCLE 1 R 2 ajustado Número parámetros BUCLE 1 R 2 ajustado Número parámetros ensayo4.b1 0.5276 24 ensayo4.poly.b1 0.5522 41 ensayo4.1.b1 0.5442 35 ensayo4.1.poly.b1 0.5586 47 BUCLE 2 BUCLE 2 ensayo4.b2 0.5227 16 ensayo4.poly.b2 0.5511 39 ensayo4.1.b2 0.5404 22 ensayo4.1.poly.b2 0.5578 34 Tabla 6. Modelos resultantes de introducir interacciones mediante bucles 1 y 2, a modelos candidatos seleccionados, de tipo lineal y polinómico. En resumen, las interacciones que parecen resultar más significativas son aquellas formadas por las covariables de predicción de la velocidad del viento y las que expresan el ciclo diario existente en el régimen de brisas; es decir, que la intensidad del viento dependerá en general, de la hora del día en que nos encontremos. El uso de bucles optimiza el tiempo invertido en la construcción de modelos, permitiendo estudiar una cantidad de interacciones entre covariables que de otro modo resultaría inviable. En particular, el bucle 2 es mucho más sencillo, tanto en eficiencia de los cálculos como en “limpieza” posterior del modelo, y obtiene modelos con menor número de parámetros y similar bondad de ajuste. 4.2.3 SELECCIÓN DE MODELOS CANDIDATOS En la Tabla 7 se han recopilado los modelos más interesantes desarrollados hasta el momento, con el fin de estudiarlos con más profundidad y seleccionar aquellos con los que seguir trabajando. Los criterios de selección utilizados han sido los expuestos en el apartado 3.2 de la metodología. De entre los modelos de tipo lineal, finalmente se han seleccionado los ensayos 4.1.b1 y 4.1.b2, que muestran los mejores resultados en términos de bondad de ajuste, desviación típica residual y desviación típica de validación cruzada (ver modelos en negrita en Tabla 7). Para guiar la elección se han realizado distintos test ANOVA, en primer lugar para la colección de ensayos lineales desarrollados a partir del ensayo4 (ver Figura A.17) y del ensayo4.1 (ver Figura A.18), para posteriormente comparar los “mejores” ensayos obtenidos en ambos casos; el test ANOVA realizado entre ellos ( ensayos 4.b2, 4.1.b1 y 4.1.b2), muestra como preferible al último (ver Figura A.19). De entre los modelos de tipo polinómico, los mejores resultados corresponden a los ensayos 4.1.poly.b1 y 4.1.poly.b2 (ver Tabla 7). En este caso, el test ANOVA muestra como más significativo al más sencillo, el ensayo4.1.poly.b2 (ver Figura A.20). 46 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín Los ensayos finalmente seleccionados se muestran en negrita en la Tabla 7. MODELO Ŝ R método CV aleatorio Ŝ R método CV consecutivo R 2 ajustado Número parámetros Ŝ R Lineal ensayo4 0.747 0.747 0.471 10 0.745 ensayo4.1 0.746 0.746 0.474 14 0.743 ensayo4.b 0.730 0.729 0.496 11 0.728 ensayo4.1.b 0.727 0.727 0.500 15 0.724 ensayo4.b1 0.708 0.708 0.528 24 0.704 ensayo4.1.b1 (mc.l1) 0.697 0.697 0.544 35 0.692 ensayo4.b2 0.711 0.711 0.523 16 0.708 ensayo4.1.b2 (mc.l2) 0.698 0.698 0.540 22 0.695 Polinómico ensayo4.poly 0.738 0.738 0.485 16 0.736 ensayo4.1.poly 0.737 0.737 0.486 16 0.734 ensayo4.poly.b 0.721 0.722 0.511 27 0.717 ensayo4.1.poly.b 0.720 0.720 0.511 22 0.716 ensayo4.poly.b1 0.692 0.691 0.552 41 0.686 ensayo4.1.poly.b1 (mc.p3) 0.688 0.687 0.559 47 0.681 ensayo4.poly.b2 0.692 0.692 0.551 39 0.686 ensayo4.1.poly.b2(mc.p4) 0.686 0.686 0.558 34 0.681 Tabla 7. Resultados de los criterios de selección estudiados para los modelos seleccionados hasta el momento: bondad de ajuste (R 2 ajustado), número de parámetros, desviación típica residual (Ŝ R ) y desviación típica resultante del proceso de validación cruzada (CV Ŝ R ) con procedimiento de ventanas aleatorias y consecutivas. En negrita se muestran los modelos seleccionados; entre paréntesis se indica el nuevo nombre asignado. Las casillas coloreadas corresponden a los mejores resultados para cada criterio de selección en estudio. El R 2 ajustado de los ensayos seleccionados oscila entre 0.559 (ensayo4.1.poly.b1, de 47 parámetros) y 0.540 (ensayo4.1.b2, de 22 parámetros). Todos ellos proceden del ensayo4.1, y a todos se han introducido las interacciones de forma automática. Se ha modificado el nombre a los modelos seleccionados, para hacer su manejo más fácil e intuitivo (ver Tabla 8). ensayo4.1.b1 mc.l1 ensayo4.1.b2 mc.l2 ensayo4.1.poly.b1 mc.p3 ensayo4.1.poly.b2 mc.p4 Tabla 8. Nuevos nombres asignados a los modelos candidatos seleccionados, donde “mc” hace referencia a modelo candidato, “l” a lineal y “p” a polinómico. En resumen, en la Tabla 7 se han coloreado las casillas que muestran los mejores resultados para cada uno de los criterios estudiados. Atendiendo a la bondad de ajuste y a la desviación típica residual, el mejor modelo candidato es mc.p3, seguido de mc.p4 (ambos polinómicos), con igual desviación residual, y que ofrece los mejores resultados en términos de desviación típica de validación cruzada. Por último, mc.l2 es interesante por su sencillez, ya que solo tiene 22 47 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín parámetros. Comparado con los otros tres, el modelo mc.l1 no parece aportar mejoras en ningún criterio.  Estudio de la composición de los modelos candidatos Los modelos candidatos mc.p3 y mc.p4 incluyen interacciones triples en su composición, muy complicadas de interpretar; también aparecen algunas interacciones “no deseadas”, por ejemplo entre términos de un armónico, que indican que los armónicos introducidos son insuficientes para recoger el comportamiento estacional y diario del régimen de vientos. Por ello, se ha estudiado la sustitución de dichas interacciones “no deseadas” por los segundos armónicos de los ciclos diarios y estacionales. El segundo armónico del ciclo diario se ha calculado de la misma forma que se calculó el primer armónico, permitiendo reproducir el comportamiento de dos ciclos durante el periodo de 24 horas. .ℎ2=ℎ 24 2∗ 2 ; .ℎ2=ℎ 24 2∗ 2 Por su parte, el segundo armónico para el ciclo estacional se ha calculado: .2=.ñ 366 2 ∗2 ; .2=.ñ 366 2 ∗2 Mediante sucesivas pruebas (no incluidas en la memoria por la extensión que éstas supondrían), se ha comprobado que a la hora de explicar el comportamiento del viento a lo largo del día, resultan necesarios también los terceros y cuartos armónicos del ciclo diario. Dichos armónicos se han calculado de la misma forma. Para los modelos candidatos lineales (mc.l1 y mc.l2) sin interacciones “no deseadas”, se ha estudiado el efecto de añadir los nuevos armónicos a los modelos, actualizando el modelo resultante; estos modelos se han renombrado a partir del modelo de origen (mc.l1.2, mc.l2.2), y se incluyen en la Tabla 9. Gracias a la introducción de los nuevos armónicos se ha conseguido eliminar ciertas interacciones triples de mc.l1, mientras el número de parámetros, bien ha aumentado en 1 (para mc.l2), bien ha disminuido (para mc.l1); por su parte, la bondad de ajuste ha mejorado ligeramente en ambos casos. También se observa cómo el porcentaje de parámetros muy significativos para el modelo aumenta en ambos casos. La introducción de los nuevos armónicos en los modelos candidatos polinómicos (mc.p3.2 y mc.p4.2) resulta beneficiosa para el modelo, pues elimina la necesidad de interacciones “no deseadas”. Además, en ambos casos permite disminuir el número de parámetros a la vez que aumenta su significación, con una bondad de ajuste y una desviación típica residual muy similar (en torno al 0.57-0.58, y entre 0.665 y 0.668 respectivamente), mientras el número de interacciones triples disminuye. También se obtienen mejores resultados en términos de desviación típica residual de validación cruzada para mc.p3 (ver casillas azules en Tabla 9). 48 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín MODELO N.par R 2 ajustado Ŝ R N.int. triples Ŝ R CV método aleatorio Ŝ R CV método consecutivo %.par.sig Lineal mc.l1 35 60% 0.544 0.692 9 0.697 0.697 mc.l1.2 27 88.9% 0.570 0.672 4 0.676 0.676 mc.l2 22 86.4% 0.540 0.695 1 0.698 0.698 mc.l2.2 23 100% 0.569 0.674 0 0.676 0.676 Polinómico mc.p3 47 61.7% 0.559 0.681 7 0.688 0.687 mc.p3.2 40 62.5% 0.579 0.665 5 0.670 0.670 mc.p4 34 85.3% 0.558 0.681 3 0.686 0.686 mc.p4.2 32 78.1% 0.575 0.668 2 0.672 0.672 Tabla 9. Estudio de la variación en los criterios de selección de los modelos candidatos, tras la introducción de nuevos armónicos de ciclo y estacionalidad. Las casillas coloreadas corresponden a los mejores resultados obtenidos para cada criterio. %.par.sig= porcentaje de parámetros con alta significación a la hora de predecir la respuesta (p-valor <0.01). N.int= número de interacciones. Los modelos en negrita corresponden a los finalmente seleccionados. Con el fin de seleccionar el modelo lineal óptimo de entre los estudiados en la Tabla 9, se ha realizado un test ANOVA que muestra como más significativo al modelo mc.l2.2 (ver Figura A. 21). Para el grupo de los modelos polinómicos, el test ANOVA no aporta gran cantidad de información, ya que éstos no están anidados (ver Figura A. 22). Antes de la selección final de los modelos óptimos, se ha estudiado el ajuste de los modelos candidatos para cada mes del periodo en estudio. Los resultados de dicho estudio vienen recogidos en la Tabla 10, donde se observa que los modelos se ajustan de forma similar a cada mes por separado, con R 2 ajustados que oscilan entre 0.63 para septiembre y 0.51 para julio. Por tanto, se considera que no será necesario desarrollar un nuevo modelo específico para ningún mes en particular. MODELO N.par R 2 ajustado R 2 ajustado Lineal Mayo Junio Julio Agosto Sept. mc.l1 35 0.544 0.519 0.537 0.514 0.558 0.603 mc.l1.2 27 0.5701 0.546 0.553 0.528 0.596 0.630 mc.l2 22 0.540 0.517 0.523 0.517 0.552 0.596 mc.l2.2 23 0.5685 0.543 0.543 0.533 0.594 0.630 Polinómico mc.p3 47 0.559 0.543 0.548 0.534 0.572 0.613 mc.p3.2 40 0.579 0.571 0.563 0.545 0.602 0.624 mc.p4 34 0.558 0.536 0.55 0.529 0.568 0.616 mc.p4.2 32 0.5754 0.562 0.561 0.540 0.595 0.623 Tabla 10. Ajuste de los modelos candidatos (con y sin nuevos armónicos) a meses por separado. Las casillas coloreadas muestran los mejores ajustes. 49 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín De entre los modelos candidatos lineales, no hay ninguno que se adapte mejor a todos los meses. Por ello, se ha seleccionado el modelo mc.l2.2, conforme al test ANOVA realizado anteriormente (ver Figura A. 21) y dada su sencillez: está formado por únicamente 23 parámetros y no contiene interacciones triples. En vistas a facilitar su análisis, se le ha denominado ms1 (“modelo seleccionado 1”). De entre los modelos candidatos polinómicos, el que mejor se ajusta, tanto al periodo en general como a los meses por separado, es mc.p3.2. Sin embargo, las diferencias con mc.p4.2 en términos de bondad de ajuste (tanto global como por meses) son mínimas. Por ello, y por ser un modelo mucho más sencillo, se ha preferido mc.p4.2, el cual se ha renombrado ms2.  Estudio comparativo de la composición de los modelos seleccionados Comparativamente, los dos modelos seleccionados utilizan 13 covariables, de las cuales 7 coinciden (ver casillas coloreadas en Tabla 11); ambos modelos incluyen covariables de alta significación. Estas covariables son de tipo predicción de la dirección del viento (cosd4), ciclo diario (sen.hora, sen.hora2, cos.hora2, sen.hora4, sen.hora4, cos.hora4) y estacionalidad (sen.dia2). Las principales diferencias entre los modelos seleccionados, son la presencia en ms2 de una covariable polinómica de segundo grado de predicción de la velocidad del viento. covariables ms1 sig. ms2 sig. intercept *** *** pred.vel 2 x4 *** 3 x2 *** x.p2 *** poly(x,2)1 *** - poly(x,2)2 *** pred.dir 2 cosd4 ** 1 cosd4 *** send3 ** - ciclo 6 sen.hora *** 7 sen.hora *** - cos.hora *** sen.hora2 *** sen.hora2 *** cos.hora2 *** cos.hora2 *** sen.hora3 ** sen.hora3 *** sen.hora4 *** sen.hora4 *** cos.hora4 ** cos.hora4 * estacionalidad 3 cos.dia *** 2 sen.dia *** sen.dia2 ** sen.dia2 *** cos.dia2 ** - TOTAL 13 13 Tabla 11. Resumen de las covariables contenidas en los modelos seleccionados ms1 y ms2, y su significación (sig.) en el modelo, a la hora de predecir la respuesta: “***” indica p-valor<0.001, “**” pvalor<0.1 y “*” p-valor <0.5. En la Tabla 12 se ha recogido el estudio de las interacciones que contienen ambos modelos: ms1 muestra únicamente 13, sin que ninguna de ellas sea triple; contiene 18, de las cuales 2 son triples. Solo 3 interacciones son comunes a ambos modelos, construidas entre predicción de la dirección y ciclo diario (send3:sen.hora, cosd4:sen.hora) y predicción de la dirección y estacionalidad (cosd4:cos.dia). De nuevo la mayor diferencia entre ambos modelos es la presencia de términos que incluyen covariables polinómicas de predicción de la velocidad del viento (de segundo y tercer grado) en ms2. 50 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín interacciones ms1 sig. ms2 sig. p.v*p.v 1 x:x2 *** 1 x3:x.p3 *** p.v*p.dir.v 2 x:send.p1 *** 3 poly(x,3)1:send3 ** x:cosd.p4 ** poly(x,3)2:send3 *** - poly(x,3)3:send3 * p.v*ciclo 2 x2:cos.hora *** 3 poly(x,3)1:sen.hora - x.p2:cos.hora *** poly(x,3)2:sen.hora - - poly(x,3)3:sen.hora ** p.v*estacionalidad 0 - 3 poly(x,3)1:cos.dia - - poly(x,3)2:cos.dia *** - poly(x,3)3:cos.dia *** p.dir.v*ciclo 3 send3:sen.hora *** 5 send3:sen.hora *** cosd4:sen.hora *** cosd4:sen.hora *** sendv:sen.hora ** send.p1:sen.hora ** - poly(x,2)1:cos.hora - - poly(x,2)3:cos.hora * p.dir.v*estacionalidad 1 cosd4:cos.dia *** 1 cosd4:cos.dia *** ciclo*estacionalidad 1 cos.hora:sen.dia *** 0 - interacciones triples 0 2 TOTAL 10 18 Tabla 12. Resumen de las interacciones contenidas en los modelos seleccionados ms1 y ms2, y su significación (sig.) en el modelo, a la hora de predecir la respuesta: “***” indica p-valor<0.001, “**” pvalor<0.1 y “*” p-valor <0.5. 4.2.4 ADECUACIÓN DE LOS MODELOS SELECCIONADOS Siguiendo la metodología desarrollada en el apartado 3.3, se ha realizado un estudio de la adecuación de los dos modelos óptimos seleccionados, basado en el análisis de sus residuos. No se ha podido realizar una validación “real” en este caso, ya que no ha sido posible disponer de una nueva BD que pudiera usarse para tal fin. Por tanto, a la hora de validar el modelo, únicamente se dispone de los resultados del proceso de validación cruzada (desviación típica residual de validación cruzada), llevada a cabo en el apartado anterior. Para los modelos seleccionados, la diferencia obtenida entre los valores de desviación típica y los de desviación típica de validación cruzada, ha sido inferior al 0.7% (ver Tabla 9) mostrando que los modelos funcionan, no solo para la BD usada en su ajuste, si no que funcionarían también para una nueva BD. ADECUACIÓN DEL MODELO CANDIDATO MS1 Los gráficos del análisis de residuos del modelo ms1 (ver Figura A.23) muestran lo siguiente: - Se verifica la hipótesis de linealidad, pues no se observa ningún patrón en el gráfico residuals vs fitted. - Se verifica la hipótesis de homocedasticidad, pues no existe evidencia de incremento o cambio de la variabilidad de la varianza. - No existen puntos influyentes que deban eliminarse del modelo (gráfico residuals vs leverage). 51 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín - Existen ciertas desviaciones de la normalidad (gráfico Normal QQ plot). Se ha descartado que el problema de la normalidad de los residuos se deba a una incorrecta transformación de la respuesta (ver estudio de la transformación Box Cox en Figura A.24). El histograma para los residuos de ms1, recogido en la Figura A.25, presenta un aspecto simétrico, con valores de media y mediana aproximadamente iguales y en torno al 0. El estudio de la distribución de los residuos estandarizados para cada mes mediante diagramas de cajas (ver Figura 29), muestra que la distribución es simétrica para todos los meses (mediana en torno al 0 y en el centro de la caja), y que por tanto, el modelo es capaz de recoger la variabilidad del comportamiento del viento de modo similar para cada mes. Los valores atípicos que superan el rango de ±3 (línea discontinua verde) no se agrupan en ningún mes en particular, siendo predominantemente de signo negativo. Se observan aproximadamente 20 datos atípicos de signo positivo en total (4 en mayo, 5 en junio, 4 en julio, 2 en agosto y 5 en septiembre). Estos residuos son de mayor interés que los de signo negativo, pues indican situaciones de alarma en las cuales el modelo no es capaz de predecir velocidades elevadas que se dan en realidad. Figura 29. Diagramas de cajas para estudiar la variación de los residuos estándar del modelo ms1 a lo largo del periodo considerado (mayo-septiembre). Los diagramas de cajas para estudiar la distribución de los residuos estándar de ms1 a lo largo del día, muestran que en concreto para el horario diurno (el horario de mayor interés, ver Figura 30), se observan 9 medidas de residuos estandarizados muy elevados, a las 13 horas (1 medida), las 16 horas (2 medidas), a las 18 horas (3 medidas) y a las 19 horas (3 medidas). También en este caso las cajas son simétricas, y la mediana se encuentra en torno al 0, mostrando un correcto comportamiento de ms1 para el horario diurno. -6 -4 -2 0 2 4 6 Residuos estandarizados vs. Meses, ms1 Meses Residuos estándar mayo junio julio agosto septiembre 52 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín Figura 30. Estudio de la variación de los diagramas de cajas para residuos estándar del modelo ms1 por horas para el horario diurno Para el horario nocturno, recogido en la Figura A. 26 del Anexo I, se observan residuos estándar mayores a +3 durante todo el periodo, con un total de 14 medidas no capturadas por el modelo (en contraste con las 9 del horario diurno). Para este horario, el modelo se ajusta peor a los datos: no se observa simetría en las cajas, y su mediana se desplaza ligeramente respecto al 0. El gráfico secuencial de residuos (recogido en la Figura A.27) no muestra ningún tipo de pauta o tendencia a lo largo del tiempo, luego no existen indicios contra la adecuación de ms1 en dicho gráfico; tampoco aportan indicios contra ella los gráficos de variables incluidas en el modelo versus residuos del modelo (ver Figura A. 28 izda. para p.v.viento y dcha. para p.dir.viento). El resultado de los gráficos crPlots para las covariables del modelo en estudio, recogido en la Figura A.29 , muestra que en general todas las covariables se adecúan bien, pues la línea verde, que indica el efecto ajustado, se superpone a la roja, que expresa el suavizado de la tendencia de los datos. Por último, el estudio de los gráficos avPlots de variable añadida (recogidos en la Figura A.30) muestra pendientes relevantes en el comportamiento de la covariable .ℎ y para la interacción :2, en cuyo caso la pendiente viene influenciada por un número escaso de puntos. La situación de estudio de mayor interés sería la existencia de grupos de residuos de signo positivo para algún parámetro, que podría dar pie a una investigación más exhaustiva sobre dichos grupos. Sin embargo, solo existen agrupaciones de residuos de signo negativo; por tanto, los gráficos avPlots tampoco aportan indicios sobre cómo mejorar el modelo ms1. ADECUACIÓN DEL MODELO CANDIDATO MS2 Para el caso del modelo ms2, que contiene covariables de tipo polinómico, los gráficos en estudio (ver Figura A.31) muestran que se recoge la linealidad y la homocedasticidad del modelo, que no existen puntos influyentes que deban eliminarse, y que existen ciertas desviaciones en -6 -4 -2 0 2 4 6 Residuos estandarizarizados vs. horario diurno, ms1 Horas Residuos estándar 8 9 10 11 12 13 14 15 16 17 18 19 59 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín Figura 37. Velocidad observada instantánea y probabilidad estimada de superar viento umbral, a partir de los modelos ms1 y ms2. Estudio de días seleccionados de forma aleatoria (21 de junio de 2012, superior; 18 de agosto de 2012, central; 11 de septiembre de 2012, inferior). 0 5 10 15 20 0 10 20 30 40 50 60 ms1. 21 junio 2012 hora v.viento (km/h) 0.25 0.5 0.75 1 P (y> umbral) umbral = 40 km/h umbral = 30 km/h umbral = 20 km/h 0 5 10 15 20 0 10 20 30 40 50 60 ms2. 21 junio 2012 hora v.viento (km/h) 0.25 0.5 0.75 1 P (y> umbral) umbral = 40 km/h umbral = 30 km/h umbral = 20 km/h 0 5 10 15 20 0 10 20 30 40 50 60 18 agosto 2012 hora v.viento (km/h) 0.25 0.5 0.75 1 P (y> umbral) umbral = 40 km/h umbral = 30 km/h umbral = 20 km/h 0 5 10 15 20 0 10 20 30 40 50 60 ms2. 18 agosto 2012 hora v.viento (km/h) 0.25 0.5 0.75 1 P (y> umbral) umbral = 40 km/h umbral = 30 km/h umbral = 20 km/h 0 5 10 15 20 0 10 20 30 40 50 60 ms1. 11 septiembre 2012 hora v.viento (km/h) 0.25 0.5 0.75 1 P (y> umbral) umbral = 40 km/h umbral = 30 km/h umbral = 20 km/h 0 5 10 15 20 0 10 20 30 40 50 60 ms2. 11 septiembre 2012 hora v.viento (km/h) 0.25 0.5 0.75 1 P (y> umbral) umbral = 40 km/h umbral = 30 km/h umbral = 20 km/h 60 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín ESTUDIO DEL PRONÓSTICO MEDIANTE TABLAS DE CONTINGENCIA Y TABLAS DE OBSERVACIONES POR DÍAS Siguiendo la metodología expuesta en el apartado 3.5, se han desarrollado dos tipos de tablas como parte del estudio del pronóstico de los modelos finalmente seleccionados, ms1 y ms2. En primer lugar se han realizado tablas de contingencia, para el estudio del comportamiento de los modelos de forma instantánea. En segundo lugar, y para estudiar su comportamiento “global” por días, es decir, si han sido capaces de pronosticar velocidades elevadas con P>0.5 en al menos un momento a lo largo del día, se han realizado un segundo tipo de tablas, que se han denominado “tablas de observaciones por días”. Recordemos que se ha considerado que el pronóstico de los modelos es correcto si son capaces de predecir un umbral de velocidad con una probabilidad superior a 0.5.  Estudio de tablas de contingencia Se han realizado tablas de contingencia para el estudio del comportamiento instantáneo de ms1 y ms2, y para umbrales de velocidad de 40 y 30 km/h. El estudio del umbral de 20 km/h se ha recogido en el Anexo III sobre estudios adicionales. Se ha obviado el estudio de umbrales inferiores de velocidad, ya que no resultan de interés para este trabajo. - Umbral > 40 km/h El umbral de 40 km/h de velocidad real se supera en 27 ocasiones a lo largo del periodo en estudio, de las cuales el modelo ms1 ha sido capaz de pronosticar 6 de forma instantánea, lo que supone un 22.22% del total de casos registrados; este porcentaje aumenta hasta el 33.33% para ms2 (ver recuadros verdes en Tabla 13). Para ambos modelos, se han dado falsas alarmas en el intervalo de velocidades de 30 a 40 km/h, en un 4.3% del total de observaciones de dicho intervalo para ms1, y en un 0.8% para ms2. No existe ninguna falsa alarma para el resto de intervalos de velocidad, en ninguno de los dos modelos (ver recuadros azules en Tabla 13). ms1 Intervalos velocidad (km/h) P(v>40 km/h) (0,10 ] (10,20 ] (20,30] (30,40] > 40 (0, 0.25] 1357 1529 600 101 18 (0.25,0.50] 0 0 1 10 3 (0.50,0.75] 0 0 0 4 2 (0.75,1] 0 0 0 1 4 ms2 (0, 0.25] 1357 1529 601 109 16 (0.25,0.50] 0 0 0 6 2 (0.50,0.75] 0 0 0 1 5 (0.75,1] 0 0 0 0 4 Tabla 13. Tabla de contingencia para pronósticos instantáneos de los modelos ms1 (superior) y ms2 (inferior), con umbral de velocidad > 40 km/h, y para la velocidad observada. En azul se muestran las falsas alarmas; en verde, las alarmas correctamente pronosticadas. 61 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín - Umbral v> 30 km/h Existen 143 observaciones en las que se ha superado el umbral de velocidad de 30 km/h durante el periodo en estudio. El modelo ms1 pronostica correctamente un 28.67% de ellas y el modelo ms2 un 27.27% (ver recuadros verdes en Tabla 14); ms1 pronostica una mayor cantidad de observaciones con P>0.75, lo cual se traduce un pronóstico más acertado. Las falsas alarmas (recuadros azules en las tablas) han sido de 19 para ms1 y de 15 para ms2. De las 27 observaciones superiores a 40 km/h, los modelos han pronosticado instantáneamente velocidades superiores a 30 km/h (por tanto, velocidades elevadas), en un 44.44% de los casos para ms1, y en un 48.15% de los casos para ms2. Es decir, este tipo de alarma puede ayudar a identificar la posible ocurrencia de “días parada”. ms1 Intervalos velocidad (km/h) P(v>30 km/h) (0,10 ] (10,20 ] (20,30] (30,40] 40 (0, 0.25] 1356 1520 538 70 9 (0.25,0.50] 0 7 46 17 6 (0.50,0.75] 1 2 16 17 3 (0.75,1] 0 0 1 12 9 ms2 (0, 0.25] 1355 1517 530 68 8 (0.25,0.50] 2 11 57 22 6 (0.50,0.75] 0 1 14 21 2 (0.75,1] 0 0 0 5 11 Tabla 14. Como Tabla 13, para probabilidad de superar los 30 km/h.  Estudio de tablas de observaciones por días En las tablas propuestas a continuación, se recogen el número de días del periodo en estudio en que se han registrado un número determinado de velocidades que superan el umbral de estudio en cada caso (40, 30 km/h), divididas entre 0, 1-2, o más de 3. Cada uno de los días en estudio se coloca en la tabla según la probabilidad con la que los modelos han pronosticado observaciones superiores al umbral en estudio, en al menos una hora a lo largo del día. De nuevo, en las tablas se han dibujado recuadros verdes, que contienen las alarmas correctamente pronosticadas por los modelos, así como recuadros azules cuando el modelo ha originado falsas alarmas. - Umbral v>40 km/h En la Tabla 15 se ha estudiado la probabilidad de superar el umbral de 40 km/h de los modelos ms1 y ms2. Existe una falsa alarma para ms1 (es decir, observaciones que el modelo ha pronosticado, serán superiores a 40 km/h, cuando en realidad no lo son), y ninguna para ms2. Para días en los que se han registrado una o dos observaciones superiores a 40 km/h (situación que ocurre en 14 ocasiones) ninguno de los modelos ha sido capaz de recoger dicho comportamiento con P>0.5. Sin embargo, ambos modelos son capaces de predecir correctamente un 66.66% de los días que han registrado entre 3 y 5 observaciones problemáticas. 62 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín ms1 Nº horas/día con observaciones > 40 km/h Intervalos P(v>40 km/h) 0 1-2 3-5 (0, 0.25] 134 11 0 (0.25,0.50] 1 3 1 (0.50,0.75] 1 0 1 (0.75,1] 0 0 1 ms2 (0, 0.25] 136 13 0 (0.25,0.50] 0 1 1 (0.50,0.75] 0 0 1 (0.75,1] 0 0 1 Tabla 15. Clasificación según probabilidad de superar los 40 km/h, de los modelos ms1 y ms2, para número de días con un número determinado de observaciones > 40 km/h Dados los resultados del estudio gráfico del funcionamiento de ms1 y ms2, los días con una o dos observaciones corresponderán en general a situaciones de tormenta. También atendiendo dichos resultados, se puede suponer que algunas de las observaciones que han superado el umbral de los 40 km/h, habrán sido pronosticadas por los modelos con al menos, P(v>30 km/h)>0.50. Se ha estudiado dicha hipótesis en la Tabla 16, donde se observa que el porcentaje de pronóstico para 1-2 observaciones aumenta, pasando de un 0 a un 28.6% para ambos modelos, mientras que para más de 3 observaciones, el porcentaje de acierto pasa del 66.66% al 100%, también para ambos modelos. En resumen, los modelos han sido capaces de pronosticar velocidades elevadas (v>30 km/h) en al menos una hora, para todos los días en que se dan 3 observaciones o más superiores al umbral de 40 km/h. En la Tabla 16, también se observan 8 falsas alarmas para ms1 y 7 para ms2; en 5 se supera realmente el umbral de los 30 km/h, aunque no para las otras 3 (ver Tabla 17). ms1 Nº horas/día con observaciones >40 km/h Intervalos P(v>30 km/h) 0 1-2 3-5 (0, 0.25] 113 5 0 (0.25,0.50] 13 5 0 (0.50,0.75] 8 3 0 (0.75,1] 0 1 3 ms2 (0, 0.25] 111 5 0 (0.25,0.50] 18 5 0 (0.50,0.75] 6 3 0 (0.75,1] 1 1 3 Tabla 16. Pronóstico de los modelos ms1 y ms2 para número de días con un número determinado de observaciones > 40 km/h, estudiando P(v>30 km/h). - Umbral v>30 km/h En la Tabla 17 se ha recogido el estudio del pronóstico de los modelos ms1 y ms2, por días y número de alarmas, para un umbral de velocidad de 30 km/h. Como se ha comentado anteriormente, ambos modelos pronostican 3 falsas alarmas para esta situación. Los modelos 63 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín son capaces de pronosticar con P>0.5, días con 1-2 observaciones con un 11.8% ms1, y con un 5.9% ms2. Para días con más de 3 observaciones, los porcentajes aumentan hasta un 62.5% en ms1 y un 56.25% en ms2. Por tanto, se ve que para este umbral de velocidad, el pronóstico de ms1 es ligeramente mejor. ms1 Nº horas/día con observaciones > 30 km/h Intervalos P(v>30 km/h) 0 1-2 3-13 (0, 0.25] 92 25 1 (0.25,0.50] 8 5 5 (0.50,0.75] 3 3 5 (0.75,1] 0 1 5 ms2 (0, 0.25] 89 25 1 (0.25,0.50] 11 7 6 (0.50,0.75] 3 2 4 (0.75,1] 0 0 5 Tabla 17. Pronóstico de los modelos ms1 y ms2 para número de días con un número determinado de observaciones > 30 km/h Como resumen final del estudio del pronóstico de los modelos óptimos seleccionados (ms1 y ms2), y en especial atendiendo a las observaciones que han superado el umbral de los 40 km/h a lo largo del periodo en estudio, se puede concluir lo siguiente: - Cuando el modelo pronostica velocidades superiores a 40 km/h, en forma de un marcado ciclo diario, puede suponerse que dicha predicción se cumplirá. - Se ha hallado que para velocidades observadas superiores a 40 km/h, los modelos han sido capaces de pronosticar, con P(v>30 km/h)>0.5, entre un 44.44% (ms1) y un 48.15% (ms2) de ellas, de forma instantánea; y entre un 28.6% (para 1-2 observaciones) y un 100% (para 3 o más), de forma diaria, es decir, para al menos una hora del día en que se han producido dichas observaciones. Por tanto, se concluye que si los modelos pronostican velocidades superiores a 30 km/h, habrá que prestar atención a dicho día y a ser posible, actualizar el modelo finalmente seleccionado con nuevas predicciones de AEMET (disponibles por ejemplo a las 0 horas del propio día pronosticado), para estudiar si éstas podrían llegar a superar el umbral de los 40 km/h. - Las dificultades que han presentado los modelos en la predicción de 1-2 observaciones diarias superiores a 40 km/h (e incluso 30 km/h), hacen pensar que éstas corresponden a eventos súbitos, tipo tormentas estivales, de rápida aparición y que producen grandes ráfagas de viento. Por su carácter local, y teniendo en cuenta que los modelos se han construido a partir del modelo de predicciones de AEMET, el cual no tiene resolución espacial y temporal suficiente, dichas situaciones de viento no podrán ser reflejadas por los modelos construidos. Por ello, se deberá prestar especial atención a los pronósticos de tormentas que operativamente pueda establecer AEMET. 64 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín 4.3 SELECCIÓN DEL MODELO ÓPTIMO, ESTUDIO DE SU CAPACIDAD OPERATIVA E IMPLEMENTACIÓN DEL PROCEDIMIENTO No se han encontrado diferencias relevantes entre los pronósticos aportados por ms1 y ms2, a excepción de que en general, ms1 calcula probabilidades más elevadas, y obtiene mejores resultados en el pronóstico de velocidades superiores a 30 km/h. Por último, ms1 funciona mejor para el día 30 de agosto (uno de los días conflictivos, ver Figura A. 40 en el Anexo I), siendo de entre los “días parada” estudiados, el que ofrece una diferencia más relevante entre ms1 y ms2. Se concluye así que la mayor complejidad del modelo ms2 (respecto a número de términos y existencia de polinomios) no ayuda a mejorar la predicción de la respuesta. Por tanto, se ha considerado a ms1 como modelo óptimo final. Para finalizar el trabajo, y conforme al apartado 3.6, se ha realizado un estudio operativo del modelo óptimo seleccionado, consistente en un estudio gráfico, y en la realización de una función en R que permitirá aplicar el modelo seleccionado operativamente a partir de las predicciones de AEMET. Para la comparación gráfica, se han realizado dos gráficos de dispersión. En la Figura 38 (izda.), se recoge el gráfico realizado entre los datos reales de velocidad (transformados en raíz) y los datos ajustados por el MRLS desarrollado a partir del modelo de AEMET, para la respuesta transformada. En la Figura 38 (dcha.) se recoge el gráfico de dispersión entre datos reales (también transformados) y datos ajustados por el modelo seleccionado (ms1). Figura 38. Gráfico de dispersión realizado entre datos iniciales y finales. A la izda. se muestra el gráfico de dispersión entre datos reales y ajustados por MRLS de AEMET; a la dcha., el correspondiente a datos reales y ajustados por ms1. Se puede verificar que el nuevo modelo de predicciones ms1 obtiene un ajuste mucho más apropiado, concentrado alrededor de la recta de identidad. Por el contrario, el ajuste del MRLS desarrollado únicamente con predicciones de AEMET resulta en un gráfico mucho más disperso 0 2 4 6 8 0 2 4 6 8 Gráfico de dispersión, v.viento vs p.v.viento AEMET v.viento predicha modelo AEMET raiz (km/h) v.viento real raiz (km/h) 02468 0 2 4 6 8 Gráfico dispersión, v.viento real vs v.viento ajustada por modelo valores v.viento ajustados modelo (raíz km/h) v.viento real (raíz km /h) 65 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín cuyos puntos no siguen la identidad, reflejando así la deficiente reproducción de la realidad observada y la escasa variabilidad que el modelo AEMET es capaz de recoger. Como finalización de este trabajo, y para permitir el uso real del modelo seleccionado ms1 se ha desarrollado un procedimiento en R. Dicho procedimiento necesita de una hoja Excel en formato .csv, en la cual se deben incluir las nuevas predicciones de AEMET a partir de las cuales se quiere obtener el pronóstico a 24 horas. En esta hoja Excel deben incluirse tres columnas: - Fecha y hora. - Predicción de AEMET de velocidad de viento para la fecha y hora determinadas. - Predicción de AEMET de dirección de viento para la fecha y hora determinadas. Para que el procedimiento funcione correctamente, las predicciones deberán corresponder a un periodo de 23 horas, desde las 3 am del día para el que se tienen las predicciones. Por supuesto, el procedimiento puede modificarse para trabajar con periodos distintos de tiempo. Como resultado, se obtendrán predicciones para el periodo entre las 7 y las 21 horas del día de interés. Los datos ajustados por el modelo (es decir, las predicciones instantáneas), vendrán incluidas en un vector. El procedimiento creará también un gráfico de pronóstico, del tipo de los estudiados en el apartado anterior (ver Figura 33). En dicho gráfico vendrán incluidas las predicciones instantáneas ajustadas por el modelo para cada hora (en el periodo entre las 7 y las 21 horas), así como las curvas de representación de probabilidad de superación de umbrales estudiados. En el Anexo III se recoge el script de R que permite el funcionamiento del procedimiento desarrollado. 66 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín 5. CONCLUSIONES Se ha realizado a lo largo de este TFG, en primer lugar, un análisis exploratorio del régimen de vientos en el vertedero de Bailín, para el periodo de mayo a septiembre de 2012, así como un estudio que demuestra la falta de adecuación a dicho régimen de vientos del modelo de predicciones de AEMET existente. Posteriormente, y conforme al objetivo general de este TFG, se han desarrollado modelos estadísticos de predicción de la velocidad del viento en forma de modelos de regresión (MR), que han permitido el downscaling del modelo de predicciones de AEMET, adaptándolo a las condiciones específicas del área y el periodo de estudio. De los análisis exploratorios se han extraído las siguientes conclusiones: 1) Se ha caracterizado en los datos observados de velocidad de viento, un ciclo diario de brisas y una estacionalidad, aunque no una tendencia. 2) Se ha establecido que el modelo de predicciones de AEMET existente funciona mejor para el periodo de invierno. Sin embargo, se ha comprobado que no resulta útil de forma directa para verano, ya que además de proporcionar valores con menor dispersión, no es capaz de reflejar la estacionalidad ni el ciclo diario existentes. Las obras de desmantelamiento y transferencia del vertedero de Bailín pueden ampliarse hasta el mes de octubre. Sin embargo, se ha decidido prescindir de dicho mes en el ajuste del nuevo modelo, ya que se considera que las predicciones de AEMET para dicho mes son apropiadas; además, los datos de octubre podrían influir en el modelo construido, limitando su capacidad predictiva para el resto de meses. 3) Los análisis exploratorios han descrito las características del viento observado y la reproducción que de ellas hacen las predicciones de AEMET. Esta fase ha proporcionado los argumentos para la construcción de la base de covariables utilizada en el desarrollo del MR óptimo buscado, y que contiene covariables de predicción de velocidad y dirección del viento (instantáneas y decaladas), así como covariables que reflejan ciclos diarios y estacionalidad, a partir de las cuales se han construido interacciones. De la metodología desarrollada para la construcción de MR, se han extraído las siguientes conclusiones: 1) El uso general de procedimientos desarrollados en R, y en particular de algoritmos para la construcción de modelos de predicción candidatos, ha optimizado el tiempo invertido, permitiendo estudiar una gran cantidad de modelos con diferentes estructuras que de otra forma habría resultado imposible. 2) Para la selección de modelos candidatos, se ha intentado hallar un equilibrio entre la complejidad de los mismos y otros criterios de selección estudiados, como la bondad de 67 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín ajuste, la desviación típica residual, la desviación típica de validación cruzada, o el resultado de los test ANOVA para la selección entre modelos anidados. Las conclusiones sobre los resultados obtenidos en la construcción de los modelos de predicción de la velocidad del viento candidatos, han sido las siguientes: 1) Las covariables que han resultado más significativas para explicar la variabilidad de la velocidad del viento real, han sido, como era de esperar por los análisis exploratorios, las covariables de predicción de la dirección del viento, tanto instantánea (extraída del modelo de AEMET) como decaladas, así como los armónicos de ciclo diario, que han aumentado la explicación de la variabilidad en un 18% y un 13% respectivamente. 2) Se ha profundizado en el estudio del comportamiento de los modelos finalmente seleccionados (ms1 lineal y ms2 polinómico), para situaciones en que se supera el umbral de los 40 km/h, velocidad a la que según el protocolo de actuación, deberían detenerse las obras. Se ha concluido que si el modelo pronostica velocidades superiores a 40 km/h, dicha predicción se cumplirá. En caso de pronosticarse velocidades superiores a 30 km/h, se deberá prestar atención a dicho día y a ser posible, actualizar el modelo finalmente seleccionado con nuevas predicciones de AEMET (disponibles por ejemplo a las 0 horas del propio día pronosticado), pues dichas observaciones son susceptibles de superar el umbral de los 40 km/h. El comportamiento de los modelos es menos satisfactorio para situaciones puntuales, como tormentas estivales, en las cuales se levantan grandes ráfagas de viento de forma súbita y de escasa duración. Por tanto, se concluye también que paralelamente al estudio del pronóstico de los modelos, se debe prestar especial atención a los pronósticos de tormentas que operativamente pueda establecer AEMET. 3) Se ha comprobado que las diferencias entre los pronósticos ofrecidos por ambos modelos estadísticos son mínimas. Por tanto, se ha concluido que los términos polinómicos, así como en general, la mayor complejidad del modelo ms2, no se traduce en una mejor predicción. 4) El modelo ms1 (lineal) finalmente seleccionado, está compuesto de 23 parámetros (13 covariables y 10 interacciones), y es capaz de explicar un 57% de la variabilidad de la velocidad del viento para el área del vertedero de Bailín. Este porcentaje representa una mejora muy importante si se tiene en cuenta que el MRLS construido a partir de las predicciones de AEMET, únicamente era capaz de explicar el 14% de dicha variabilidad. A la hora de juzgar este porcentaje, se debe tener en cuenta también la dificultad que supone trabajar con datos meteorológicos, cuyas características hacen de su predicción una ardua tarea. 5) Para que el modelo construido pueda ser utilizado, se ha desarrollado un procedimiento en R, que mediante la introducción de las predicciones de AEMET de interés (para un periodo de 23 horas) en una hoja Excel, es capaz de ofrecer un estudio gráfico del pronóstico del modelo, en términos de predicciones instantáneas y de probabilidad de ocurrencia de distintos umbrales de velocidad del viento. 68 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín Para resumir, a lo largo de este TFG se ha desarrollado una metodología de construcción de un modelo de predicción de la velocidad del viento, que ha permitido el downscaling de la predicción AEMET, el cual se ha aplicado en Bailín con resultados satisfactorios. El modelo desarrollado pretende ser una herramienta más en la lucha contra la contaminación de áreas adyacentes al vertedero de Bailín a causa de la dispersión atmosférica, lo cual incrementaría la problemática ambiental existente.  Posibles trabajos futuros Teniendo en cuenta el carácter herramental que poseen los modelos de regresión para distintas disciplinas, la metodología desarrollada puede servir para resolver otro tipo de problemas, en los cuales se necesite predecir cualquier respuesta a partir de un conjunto de covariables predictoras. Para este TFG en concreto, se han resumido una serie de recomendaciones en vistas a una posible realización de trabajos futuros que sigan la línea del estudio realizado. Las recomendaciones que se sugieren son las siguientes: - En estudios incluidos en el Anexo III, se ha visto cómo la temperatura real muestra una clara relación con la velocidad del viento. Lamentablemente, a la hora de construir el modelo no se ha podido disponer de datos de predicción de la temperatura. La introducción de esas (u otras) variables atmosféricas para la zona, permitirían ajustar mejor el modelo, obteniendo pronósticos más acertados. - Se cree que trabajando con predicciones AEMET actualizadas a más corto plazo, se obtendrían predicciones más acordes a la realidad (en especial en lo que a velocidades elevadas se refiere), tanto de forma instantánea como para el estudio del pronóstico diario. Es decir, si el modelo trabajase operativamente, nutriéndose de forma continua de nuevas predicciones de AEMET actualizadas, podría mejorar las probabilidades obtenidas. 75 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín Figura A. 7. Test ANOVA entre ensayos 4 y 4.1 (en orden de aparición). El ensayo4.1 resulta altamente significativo. Figura A.8. Estudio gráfico de los grados más apropiados para las covariables de predicción de viento del ensayo4, mediante la función gam en R. 76 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín Figura A.9. Estudio de modelos anidados para ensayo4, analizando distintos grados para las covariables de predicción de la velocidad del viento. Los modelos se han creado mediante la función gam en R, y comparado con la función anova, para hallar el más significativo (marcado en amarillo). Figura A.10. Resultado del test ANOVA aplicado para estudiar el grado de x4 en ensayo4.1, obteniendo como resultado que x4 se adapta mejor al grado 1. Figura A.11 Test ANOVA realizado entre ensayo4 y ensayo4.poly. Como resultado, se obtiene como más significativo al ensayo4.poly. 77 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín Figura A.12 Test ANOVA realizado entre ensayo4.1 y ensayo4.1.poly. Como resultado, se obtiene como más significativo al ensayo4.1.poly. Figura A.13 Resultado del test ANOVA entre los 8 modelos lineales a los que se ha introducido interacciones de forma manual. Los modelos 1 a 8 siguen el orden de aparición de la columna de modelos lineales en la Tabla 5. Figura A.14. Resultado del test ANOVA entre los 8 modelos polinómicos a los que se ha introducido interacciones de forma manual. Los modelos aquí representados (del 1 al 8) siguen el orden de aparición de los modelos polinómicos en la Tabla 5. 78 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín Figura A.15 Resultados del test ANOVA realizado para los modelos lineales desarrollados mediante b1 y b2, estudiados en el orden en que aparecen en la Tabla 6. El modelo más significativo es el ensayo4.1.b2 Figura A.16 Resultados del test ANOVA realizado para los modelos polinómicos desarrollados mediante b1 y b2, estudiados en el orden en que aparecen en la Tabla 6. El modelo más significativo, de entre los que han resultado comparables por estar anidados, es el ensayo4.poly.b2. No se ha podido comparar con ensayo4.1.poly.b2. Figura A.17. Resultado del test ANOVA para ensayos construidos a partir del ensayo4 (4, 4.b, 4.b1 y 4.b2); resulta más significativo el ensayo4.b2. Figura A.18 Resultado del test ANOVA para ensayos construidos a partir del ensayo4.1; de entre los modelos anidados, resulta más significativo ensayo4.1.b2. 79 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín Figura A.19 Resultado del test ANOVA entre ensayo4.b2, ensayo4.1.b1 y ensayo4.1.b2 (por orden de aparición). Figura A.20 Resultados del test ANOVA para los ensayos 4.1.poly.b1 y 4.1.poly.b2, mostrando como más significativo al último. Figura A. 21 Resultados del test ANOVA para los modelos lineales candidatos (en orden de aparición en Tabla 9), para medir la significación de los nuevos armónicos introducidos. Figura A. 22 Resultados del test ANOVA para los modelos polinómicos candidatos (en orden de aparición en Tabla 9), para medir la significación de los nuevos armónicos introducidos. 80 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín Figura A.23. Resultados de la función plot en R, para modelo ms1. Comprobación de hipótesis de residuos. Figura A.24. Resultado del estudio de la transformación Box Cox para ms1. La transformación correcta de la respuesta es la transformación raíz cuadrada. 2 3 4 5 6 7 -4 -2 0 2 4 Fitted values Residuals Residuals vs Fitted 2859 2593 3499 -2 0 2 -4 0 2 4 6 Theoretical Quantiles Standardized residuals Normal Q-Q 2859 2593 3499 2 3 4 5 6 7 0.0 1.0 2.0 Fitted values Standardized residuals Scale-Location 2859 2593 3499 0.000 0.010 0.020 -6 -2 2 4 6 Leverage Standardized residuals Cook's distance Residuals vs Leverage 34992107 1472 -1.0 -0.5 0.0 0.5 1.0 -28000 -26000 -24000 -22000 λ log-Likelihood 95% 81 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín Figura A.25. Histograma de residuos para ms1. Figura A. 26. Diagramas de cajas para estudiar la variación de los residuos estándar del modelo ms1 a lo largo del horario nocturno (entre las 20 y las 7 horas, para cada hora por separado). Histograma ms1 Residuos ms1 Frecuencia -6 -4 -2 0 2 4 6 0 500 1000 1500 media = -1.00e-17 mediana = -0.018 -6 -4 -2 0 2 4 6 Residuos estandarizados vs. horario nocturno, ms1 Horas Residuos estándar 20 21 22 23 0 1 2 3 4 5 6 7 82 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín Figura A.27. Gráfico residuos vs tiempo para ms1. No se observa ningún patrón especial. Figura A. 28. Gráfico de residuos del modelo ms1 versus covariables de predicción de la velocidad del viento del modelo de AEMET (izda.) y de predicción de la dirección del viento del modelo AEMET (dcha.). 0123456 -6 -4 -2 0 2 4 6 Residuos ms1 vs predicción velocidad (raiz) raiz p.v.viento (km/h) residuos ms1 0 50 100 150 200 250 300 350 -6 -4 -2 0 2 4 6 Residuos ms1 vs predicción dirección p.dir.viento (grados) Residuos Ms1 83 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín Figura A.29. Gráficos crPlots para el estudio de la adecuación de las covariables del modelo ms1. Continúa. 84 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín Figura A.29. Gráficos crPlots para el estudio de la adecuación de las covariables del modelo ms1. Figura A.30. Gráficos avPlots para el modelo ms1. Continúa. 91 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín Figura A.38 Gráficos avPlots para el modelo ms2. Continúa. 92 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín Figura A.38. Gráficos avPlots para el modelo ms2. Continúa. 93 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín Figura A.38. Gráficos avPlots para el modelo ms2. ESTUDIO GRÁFICO DEL PRONÓSTICO DE LOS MODELOS SELECCIONADOS MS1 Y MS2  Estudio de “días parada” Mes Día Hora mayo 17 18 21 11 21 12 21 13 21 14 21 15 junio 11 16 12 16 18 19 julio 1 17 26 22 27 18 agosto 25 17 28 2 30 18 31 17 septiembre 1 16 1 19 3 16 12 14 12 15 12 16 13 12 13 13 13 14 13 15 23 18 Tabla A.39. Tabla resumen de “días parada”, considerando como tales, aquellos que presentan observaciones de velocidad de viento superiores a 40 km/h. 94 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín  “Días parada” en los que el modelo pronostica correctamente “Días parada” con observaciones problemáticas que el modelo pronostica correctamente “Días parada” con observaciones aisladas 0 5 10 15 20 0 10 20 30 40 50 60 ms1. 13 septiembre 2012 hora v.viento (km/h) 0.25 0.5 0.75 1 P (y> umbral) umbral = 40 km/h umbral = 30 km/h 0 5 10 15 20 0 10 20 30 40 50 60 ms2. 13 septiembre 2012 hora v.viento (km/h) 0.25 0.5 0.75 1 P (y> umbral) umbral = 40 km/h umbral = 30 km/h 0 5 10 15 20 0 10 20 30 40 50 60 ms1. 11 junio 2012 hora v.viento (km/h) 0.25 0.5 0.75 1 P (y> umbral) umbral = 40 km/h umbral = 30 km/h umbral = 20 km/h 0 5 10 15 20 0 10 20 30 40 50 60 ms2. 11 junio 2012 hora v.viento (km/h) 0.25 0.5 0.75 1 P (y> umbral) umbral = 40 km/h umbral = 30 km/h umbral = 20 km/h 95 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín 0 5 10 15 20 0 10 20 30 40 50 60 ms1. 18 junio 2012 hora v.viento (km/h) 0.25 0.5 0.75 1 P (y> umbral) umbral = 40 km/h umbral = 30 km/h umbral = 20 km/h 0 5 10 15 20 0 10 20 30 40 50 60 ms2. 18 junio 2012 hora v.viento (km/h) 0.25 0.5 0.75 1 P (y> umbral) umbral = 40 km/h umbral = 30 km/h umbral = 20 km/h 0 5 10 15 20 0 10 20 30 40 50 60 ms1. 27 julio 2012 hora v.viento (km/h) 0.25 0.5 0.75 1 P (y> umbral) umbral = 40 km/h umbral = 30 km/h umbral = 20 km/h 0 5 10 15 20 0 10 20 30 40 50 60 ms2. 27 julio 2012 hora v.viento (km/h) 0.25 0.5 0.75 1 P (y> umbral) umbral = 40 km/h umbral = 30 km/h umbral = 20 km/h 0 5 10 15 20 0 10 20 30 40 50 60 ms1. 31 agosto 2012 hora v.viento (km/h) 0.25 0.5 0.75 1 P (y> umbral) umbral = 40 km/h umbral = 30 km/h 0 5 10 15 20 0 10 20 30 40 50 60 ms2. 31 agosto 2012 hora v.viento (km/h) 0.25 0.5 0.75 1 P (y> umbral) umbral = 40 km/h umbral = 30 km/h 96 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín “Días parada” en que el modelo pronostica velocidades elevadas 0 5 10 15 20 0 10 20 30 40 50 60 ms1. 23 septiembre 2012 hora v.viento (km/h) 0.25 0.5 0.75 1 P (y> umbral) umbral = 40 km/h umbral = 30 km/h umbral = 20 km/h 0 5 10 15 20 0 10 20 30 40 50 60 ms2. 23 septiembre 2012 hora v.viento (km/h) 0.25 0.5 0.75 1 P (y> umbral) umbral = 40 km/h umbral = 30 km/h umbral = 20 km/h 0 5 10 15 20 0 10 20 30 40 50 60 ms1. 12 junio 2012 hora v.viento (km/h) 0.25 0.5 0.75 1 P (y> umbral) umbral = 40 km/h umbral = 30 km/h 0 5 10 15 20 0 10 20 30 40 50 60 ms2. 12 junio 2012 hora v.viento (km/h) 0.25 0.5 0.75 1 P (y> umbral) umbral = 40 km/h umbral = 30 km/h 0 5 10 15 20 0 10 20 30 40 50 60 ms1. 12 septiembre 2012 hora v.viento (km/h) 0.25 0.5 0.75 1 P (y> umbral) umbral = 40 km/h umbral = 30 km/h 0 5 10 15 20 0 10 20 30 40 50 60 ms2. 12 septiembre 2012 hora v.viento (km/h) 0.25 0.5 0.75 1 P (y> umbral) umbral = 40 km/h umbral = 30 km/h 97 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín  “Días parada” en los que el modelo no pronostica correctamente, sin causa aparente 0 5 10 15 20 0 10 20 30 40 50 60 ms1. 1 septiembre 2012 hora v.viento (km/h) 0.25 0.5 0.75 1 P (y> umbral) umbral = 40 km/h umbral = 30 km/h umbral = 20 km/h 0 5 10 15 20 0 10 20 30 40 50 60 ms2. 1 septiembre 2012 hora v.viento (km/h) 0.25 0.5 0.75 1 P (y> umbral) umbral = 40 km/h umbral = 30 km/h umbral = 20 km/h 0 5 10 15 20 0 10 20 30 40 50 60 ms1. 3 septiembre 2012 hora v.viento (km/h) 0.25 0.5 0.75 1 P (y> umbral) umbral = 40 km/h umbral = 30 km/h umbral = 20 km/h 0 5 10 15 20 0 10 20 30 40 50 60 ms2. 3 septiembre 2012 hora v.viento (km/h) 0.25 0.5 0.75 1 P (y> umbral) umbral = 40 km/h umbral = 30 km/h umbral = 20 km/h 98 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín  Día con diferencias más visibles entre ms1 y ms2 Figura A. 40 Estudio del pronóstico de ms1 y ms2 para el 30 de agosto de 2012, día con diferencias más significativas en el pronóstico de ambos modelos. 0 5 10 15 20 0 10 20 30 40 50 60 ms1. 30 agosto 2012 hora v.viento (km/h) 0.25 0.5 0.75 1 P (y> umbral) umbral = 40 km/h umbral = 30 km/h 0 5 10 15 20 0 10 20 30 40 50 60 ms2. 30 agosto 2012 hora v.viento (km/h) 0.25 0.5 0.75 1 P (y> umbral) umbral = 40 km/h umbral = 30 km/h umbral = 20 km/h 99 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín ANEXO II. CONCEPTOS ESTADÍSTICOS DE INTERÉS LAS ROSAS DE LOS VIENTOS La rosa de los vientos (ver ejemplo en la Figura A.41) es una herramienta que permite resumir de forma gráfica el comportamiento del viento, ya que representa la frecuencia con que se ha registrado una determinada velocidad en cada partición de la dirección del viento (en este caso definida por 12 rumbos), mediante códigos de amplitud y de color. El código de color muestra los intervalos de velocidad considerados, mientras la amplitud de la franja coloreada muestra la frecuencia a la que dichos intervalos se han dado, para cada dirección. Figura A.41 Ejemplo de rosa de vientos La realización de rosas de vientos permitirá comparar gráficamente el comportamiento del viento para las BD real y predicha por el modelo de AEMET. CONCEPTOS SOBRE CORRELACIÓN El coeficiente de correlación de Pearson indica cuál es la relación lineal entre un par de variables para el instante de tiempo t 0 , es decir, indica la correlación instantánea entre ellas. Sin embargo, en ocasiones interesa estudiar también cuál es la correlación para cada instante de tiempo del periodo en estudio, es decir, la correlación cruzada. Este análisis se realiza cuando se quiere calcular cuál sería la correlación para instantes decalados. La autocorrelación permite estudiar la relación existente entre los valores de una misma variable, observados con una diferencia de k unidades de tiempo, denominada “retardo”. El estudio de la autocorrelación de las variables aporta información sobre su comportamiento a lo largo del tiempo; es decir, permite estudiar por ejemplo, si los datos anteriores o posteriores de una variable en el tiempo, podrían ser útiles a la hora de explicar el comportamiento de dicha variable. 100 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín Así, la existencia de autocorrelación puede expresar la existencia de estacionalidad, lo cual es de esperar si, como en nuestro caso, se trabaja con series de tiempo horarias (Pepió, 2011). El gráfico que permite estudiar la correlación, tanto entre un par de variables como para una misma variable, se denomina correlograma. El estudio del correlograma resulta útil cuando, por ejemplo, se cree que puede existir una correlación retardada (es decir, una correlación mayor con otro instante de tiempo), una autocorrelación, o para estudiar la amplitud de los ciclos existentes. El correlograma se compone de:  En el eje x, el periodo considerado para nuestros datos. En nuestro caso, este periodo será de 24 horas.  En el eje y, bien la correlación cruzada en caso de estudiar dos series de datos, bien la autocorrelación en caso de estudiar una única variable. INTRODUCCIÓN A LOS MODELOS DE REGRESIÓN LINEAL Según el objetivo general de este TFG, se ha planteado el desarrollo de una serie de modelos de regresión (MR), lineales y no lineales, que permiten explicar y predecir el comportamiento de una variable cuantitativa (<), llamada respuesta, a partir de una o un conjunto de variables regresoras o predictoras (=), suponiendo que la relación entre ellas es lineal o linealizable. Además, según James et al. (2013), los MR permiten contestar a cuestiones que se plantean en el trabajo, tales como: - ¿Existe una relación entre predictores y respuesta? ¿Cuál es su forma? (es lineal, no lineal). - A la hora de predecir la respuesta, ¿resulta útil todo el conjunto de los predictores, o es solamente un subconjunto de ellos el que ayuda a explicarla? - ¿Se ajusta el modelo de regresión a nuestro conjunto de datos?, es decir, ¿se ajustan a la realidad las predicciones de nuestro modelo? ¿cómo de precisa es dicha predicción? Los modelos de regresión lineal (MRL) pueden ser simples (MRLS), o múltiples (MRLM), dependiendo del número de variables regresoras en estudio. Así, si únicamente se dispone de una variable regresora, el modelo será simple y seguirá la expresión recogida en Ec.2, llamada ecuación de regresión, donde < es la variable respuesta, = la variable predictora, > representa el residuo o error de la regresión, '  es el punto de corte de la recta con el eje de abscisas y ' ( representa la pendiente de la recta. '  y ' ( se denominan coeficientes de regresión. < = '  +' ( =+ > (Ec. 2) Sin embargo, normalmente no se dispone de una única variable regresora en estudio, sino de un conjunto de ellas. En ese caso, el modelo en estudio se denomina de regresión lineal múltiple y es una generalización del modelo anterior: < = '  +' ( = ( +'  =  +⋯+' $ = $ + > (Ec. 3) El fin del MRLM sigue siendo el mismo: usar el conjunto de variables regresoras = para prededir la respuesta <. En este caso, la línea de regresión se reemplaza por un hiperplano. 107 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín ANEXO III. ESTUDIOS ADICIONALES DESCOMPOSICIÓN DE SERIES TEMPORALES EN R La Figura A. 43 muestra el estudio de la descomposición de la serie temporal en estudio (BD.verano) en R. En ella se puede observar fácilmente la existencia de ciclo y estacionalidad (seasonal), aunque no se aprecia ningún patrón en la tendencia (trend). Por tanto, se concluye que la serie de datos reales en estudio no posee una tendencia, probablemente debido a que el periodo en estudio no es lo suficientemente largo como para apreciarla. Figura A. 43. Descomposición de la BD.max2, como serie de datos temporal. En esta descomposición se puede observar una estacionalidad, pero no una tendencia clara. ESTUDIO DE LOS DATOS DE TEMPERATURA REAL  Estudio de la relación instantánea Se ha realizado un estudio paralelo al desarrollo del trabajo, con el fin de averiguar si existe alguna relación entre los datos de velocidad real y temperatura, pues vista la estacionalidad diaria del régimen de vientos, se cree que también existirá dicha relación lineal entre ambas covariables, ya que el régimen de temperaturas también ofrece diferencias notables entre sus medidas registradas a lo largo del día. Los datos de temperatura real utilizados se han seleccionado con el mismo procedimiento usado para la selección de los datos de dirección correspondientes a las velocidades de viento de interés (ver descripción de P2 en apartado 2.1.2 de la metodología). El gráfico de dispersión entre la velocidad real y la temperatura real, recogido en la Figura A. 44, muestra la existencia de una clara relación lineal entre ambas covariables; en términos generales, conforme aumenta la velocidad del viento aumenta la temperatura, y viceversa (línea azul gráfico). Esta situación se cumple bastante bien hasta velocidades de 25-30 km/h, aunque no para 108 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín velocidades mayores (datos incluidos en el círculo verde en el gráfico), que seguramente, no sigan un comportamiento “normal”, si no que se trate de eventos puntuales. Figura A. 44 Gráfico de dispersión entre las covariables velocidad real de viento y temperatura.  Estudios de correlación La correlación instantánea entre las covariables temperatura real y velocidad del viento real (de la BD.verano) es de 0.4466, lo cual indica que existe una relación lineal importante entre ambas variables. Además, al analizar el gráfico de correlación cruzada (ver Figura A. 45), se observa que ambas variables siguen un ciclo sinusoidal, donde la tendencia de dicha curva cambia cada 12 horas. Figura A. 45 Correlograma entre covariables de temperatura y velocidad real del viento, para el periodo en estudio. -30 -20 -10 0 10 20 30 -0.4 -0.2 0.0 0.2 0.4 Lag ACF Correlación entre v.viento y temperatura 109 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín Los resultados obtenidos mediante el análisis exploratorio de los datos reales de temperatura, hacen pensar que el ciclo diario existente en la velocidad del viento, puede estar altamente influenciado por la variación diaria de la temperatura. Sin embargo, no se ha podido disponer de datos de predicción de temperatura y por tanto, no se ha podido completar con ellos el modelo estadístico de predicción desarrollado. ESTUDIO DEL PRONÓSTICO DE LOS MODELOS SELECCIONADOS. UMBRAL DE 20 KM/H Ya que el interés del trabajo radica en las velocidades de viento más elevadas (superiores a 30 km/h), se ha incluido en este apartado el estudio realizado para el pronóstico de los modelos ms1 y ms2 finalmente seleccionados, para el umbral de velocidades superiores a 20 km/h, ya que como se ha dicho, tienen un menor interés para este trabajo. Se ha llevado a cabo el mismo proceso realizado en el apartado 4.2.5. En primer lugar, se han llevado a cabo tablas de contingencia para comparar el comportamiento del modelo forma instantánea, tomando como umbral v>20 km/h. En este caso, el modelo ms1 ofrece un pronóstico instantáneo acertado en un 51.5% de los casos para ms1, para el intervalo entre 20 y 30 km/h, y en un 51.1% de los casos para ms2 (ver recuadros verdes en Tabla A. 46). ms1 Intervalos velocidad (km/h) P(v>20 km/h) (0,10 ] (10,20 ] (20,30] (30,40] 40 (0, 0.25] 1325 1089 156 19 3 (0.25,0.50] 23 321 166 13 4 (0.50,0.75] 8 102 191 31 2 (0.75,1] 1 17 88 53 18 ms1 (0, 0.25] 1318 1096 145 13 1 (0.25,0.50] 30 319 185 15 5 (0.50,0.75] 7 94 184 40 2 (0.75,1] 2 20 87 48 19 Tabla A. 46 Como Tabla 13, para probabilidad de superar los 20 km/h En segundo lugar, se han construido tablas que incluyen el número de observaciones/día, y la capacidad de los modelos para pronosticar el umbral en estudio, en al menos una ocasión a lo largo del día en que se han dado las observaciones. Para el umbral de velocidad de 20 km/h, ambos modelos pronostican falsas alarmas (1 en el caso de ms1 y 2 en el caso de ms2). Respecto al porcentaje de aciertos en días con 1-2 observaciones, éste supone un 26.47% para ambos modelos. El porcentaje de aciertos para días con más de 3 observaciones (y hasta 19), es de 78.4% para ms1 y 82% para ms2 (ver Tabla A. 47). Por tanto, para el umbral de 20 km/h, es el modelo ms2 el que funciona ligeramente mejor. 110 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín ms1 Nº horas/día con observaciones >20 km/h Intervalos P(v>20 km/h) 0 1-2 3-19 (0, 0.25] 3 4 2 (0.25,0.50] 6 19 22 (0.50,0.75] 1 8 45 (0.75,1] 0 1 42 ms2 (0, 0.25] 1 2 0 (0.25,0.50] 7 21 20 (0.50,0.75] 2 7 46 (0.75,1] 0 2 45 Tabla A. 47 Pronóstico de los modelos ms1 y ms2 para número de días con un número determinado de observaciones > 20 km/h SCRIPT DE R QUE PERMITE OBTENER PREDICCIONES A 24 HORAS, APLICANDO EL MODELO ÓPTIMO DESARROLLADO (MS1) #### IMPORTANTE: ANTES DE EMPEZAR, HAY QUE INSTALAR LOS PAQUETES lubridate, chron y #### modeest setwd('C:/') # Aquí es donde se encuentra el archivo Excel al cual introducir predicciones de # AEMET, para el día de interés (entre las 3 am y la 1 am). El archivo se debe # llamar prediccion y tener formato .csv. Los decimales de la velocidad del viento # vendrán especificados con puntos. Fecha y hora se especificarán en la misma # columna fecha.hora, predicción de la velocidad del viento vendrá en la # columna p.v.viento y predicción de la dirección del viento en la columna # p.dir.viento prevision.hoy <- read.csv("prediccion.csv", header=T, sep=';') prevision.hoy <- as.data.frame(prevision.hoy) attach(prevision.hoy) # 1. Se separa “hora” para que el programa la pueda leer auxiliar.fecha<-strptime(fecha.hora, "%d/%m/%Y %H:%M") require(lubridate) # devuelve mes como número require(chron) require(modeest) prevision.hoy$año<-year(auxiliar.fecha) prevision.hoy$mes<-month(auxiliar.fecha) prevision.hoy$dia<-day(auxiliar.fecha) prevision.hoy$hora<-hour(auxiliar.fecha) # 2. Se crea un data.frame que contenga todas las covariables necesarias que el modelo necesita # 1. Se crean covariables post y lag, de v.viento y de p.v.viento 111 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín # Para v.viento attach(prevision.hoy) auxiliar.nombre.variable <- 'p.v.viento' auxiliar.variable.interes <- prevision.hoy[[auxiliar.nombre.variable ]] for(i.lag in 1:4) { auxiliar.nombre.variable.interes <- paste(auxiliar.nombre.variable, 'lag', i.lag, sep='.', collapse='') auxiliar.variable<-c(rep(NA,i.lag), auxiliar.variable.interes[1:(length(auxiliar.variable.interes)-i.lag)]) prevision.hoy <- cbind(prevision.hoy, auxiliar.variable) names(prevision.hoy)[length(names(prevision.hoy))] <- auxiliar.nombre.variable.interes } for(i.post in 1:4) { auxiliar.nombre.variable.interes <- paste(auxiliar.nombre.variable, 'post', i.post, sep='.', collapse='') auxiliar.variable <- c(auxiliar.variable.interes[(1+i.post):length(auxiliar.variable.interes)],rep(NA, i.post)) # movemos los datos prevision.hoy <- cbind(prevision.hoy, auxiliar.variable) names(prevision.hoy)[length(names(prevision.hoy))] <- auxiliar.nombre.variable.interes } # Para p.v.viento auxiliar.nombre.variable <- 'p.dir.viento' auxiliar.variable.interes <- prevision.hoy[[auxiliar.nombre.variable ]] for(i.lag in 1:4) { auxiliar.nombre.variable.interes <- paste(auxiliar.nombre.variable, 'lag', i.lag, sep='.', collapse='') auxiliar.variable <-c(rep(NA,i.lag), auxiliar.variable.interes[1:(length(auxiliar.variable.interes)-i.lag)]) prevision.hoy <- cbind(prevision.hoy, auxiliar.variable) names(prevision.hoy)[length(names(prevision.hoy))]<- auxiliar.nombre.variable.interes } for(i.post in 1:4) { auxiliar.nombre.variable.interes <- paste(auxiliar.nombre.variable, 'post', i.post, sep='.', collapse='') auxiliar.variable <- c(auxiliar.variable.interes[(1+i.post):length(auxiliar.variable.interes)],rep(NA, i.post)) # movemos los datos 112 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín prevision.hoy <- cbind(prevision.hoy, auxiliar.variable) names(prevision.hoy)[length(names(prevision.hoy))]<- auxiliar.nombre.variable.interes } prevision.hoy # comprobamos que se han añadido las nuevas covariables # Se cambia el el nombre a las covariables anteriores, para facilitar su manejo # velocidad detach() attach(prevision.hoy) x <- p.v.viento x1 <- p.v.viento.lag.1 x2 <- p.v.viento.lag.2 x3 <- p.v.viento.lag.3 x4 <- p.v.viento.lag.4 x.p1 <- p.v.viento.post.1 x.p2 <- p.v.viento.post.2 x.p3 <- p.v.viento.post.3 x.p4 <- p.v.viento.post.4 # dirección dx <- p.dir.viento d1 <- p.dir.viento.lag.1 d2 <- p.dir.viento.lag.2 d3 <- p.dir.viento.lag.3 d4 <- p.dir.viento.lag.4 d.p1 <- p.dir.viento.post.1 d.p2 <- p.dir.viento.post.2 d.p3 <- p.dir.viento.post.3 d.p4 <- p.dir.viento.post.4 # Se crean los armónicos de dirección de viento cosdv <- cos(dx/360 * 2 * pi) sendv <- sin(dx/360 * 2 * pi) cosd1 <- cos(d1/360 * 2 * pi) send1 <- sin(d1/360 * 2 * pi) cosd2 <- cos(d2/360 * 2 * pi) send2 <- sin(d2/360 * 2 * pi) cosd3 <- cos(d3/360 * 2 * pi) send3 <- sin(d3/360 * 2 * pi) cosd4 <- cos(d4/360 * 2 * pi) send4 <- sin(d4/360 * 2 * pi) cosd.p1 <- cos(d.p1/360 * 2 * pi) send.p1 <- sin(d.p1/360 * 2 * pi) 113 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín cosd.p2 <- cos(d.p2/360 * 2 * pi) send.p2 <- sin(d.p2/360 * 2 * pi) cosd.p3 <- cos(d.p3/360 * 2 * pi) send.p3 <- sin(d.p3/360 * 2 * pi) cosd.p4 <- cos(d.p4/360 * 2 * pi) send.p4 <- sin(d.p4/360 * 2 * pi) prediccion.hoy <- cbind(mes, dia, hora, x,x1,x2,x3,x4,x.p1,x.p2,x.p3,x.p4, sendv, send1, send2, send3, send4, send.p1, send.p2, send.p3, send.p4, cosdv, cosd1, cosd2, cosd3, cosd4, cosd.p1, cosd.p2, cosd.p3, cosd.p4) summary(prediccion.hoy) # Se comprueba que la BD se ha creado correctamente prediccion.hoy <- as.data.frame(prediccion.hoy) detach() attach(prediccion.hoy) # Se crean armónicos de estacionalidad y ciclos # Estacionalidad dias.transcurridos <- 121 dia.estudio <- rep(dias.transcurridos + prediccion.hoy$dia + (prediccion.hoy$mes - 5)*31) sen.dia <- sin(dia.estudio *2*pi/365) cos.dia <- cos(dia.estudio *2*pi/365) sen.dia2 <- sin(dia.estudio *4*pi/365) cos.dia2 <- cos(dia.estudio *4*pi/365) # Ciclo sen.hora <- sin(prediccion.hoy$hora*2*pi/24) cos.hora <- cos(prediccion.hoy$hora*2*pi/24) sen.hora2 <- sin(prediccion.hoy$hora*4*pi/24) cos.hora2 <- cos(prediccion.hoy$hora*4*pi/24) sen.hora3 <- sin(prediccion.hoy$hora*8*pi/24) cos.hora3 <- cos(prediccion.hoy$hora*8*pi/24) sen.hora4 <- sin(prediccion.hoy$hora*16*pi/24) cos.hora4 <- cos(prediccion.hoy$hora*16*pi/24) prediccion.hoy <- cbind(prediccion.hoy, sen.dia, cos.dia, sen2.dia, cos2.dia, sen.hora, cos.hora, sen2.hora, cos2.hora, sen3.hora, cos3.hora, sen4.hora, cos4.hora) 114 Desarrollo de un modelo estadístico de predicción de la velocidad del viento para el área del vertedero de Bailín # 3. Se aplica el modelo a nuevo conjunto de datos, con función predict() # Se carga área de trabajo BD.final.RData (donde está el modelo óptimo ms1 # guardado). Es importante cargarla en este momento y no antes, ya que de lo # contrario la base de datos no se creará correctamente prediccion.prueba <- predict.lm(ms1, newdata=prediccion.hoy, se.fit=TRUE, type='response', terms=NULL) # 4. Se calculan los pronósticos, para umbrales de velocidad de 40, 30, 20 y 10 km/h horas.utiles <- c(7:21) auxiliar.p.40 <- pnorm(sqrt(40) , mean=na.omit(prediccion.prueba$fit), sd=prediccion.prueba$residual.scale) auxiliar.p.30 <- pnorm(sqrt(30) , mean=na.omit(prediccion.prueba$fit) , sd=prediccion.prueba$residual.scale) auxiliar.p.20 <- pnorm(sqrt(20) , mean=na.omit(prediccion.prueba$fit), sd=prediccion.prueba$residual.scale) auxiliar.p.10 <- pnorm(sqrt(10) , mean=na.omit(prediccion.prueba$fit), sd=prediccion.prueba$residual.scale) plot(horas.utiles, (na.omit(prediccion.prueba$fit))^2, ylim=c(0,60), xlab='hora',ylab='v.viento (km/h)', pch=19, main='Gráfico predicciones instantáneas y pronósticos') points(prevision.hoy$v.viento[horas.utiles], pch=19, col='blue') abline(h=c(10, 20,30,40), lty=2) lines(horas.utiles, (1-auxiliar.p.40)*60,col='red', lwd=2) abline(h=c(0.25,0.5,0.75)*60,col='red',lty=2) axis(4,at=60*c(0.25,0.5,0.75,1),labels=c(0.25,0.5,0.75,1),tick=T) mtext(side=4,text='P (y> umbral)') lines(horas.utiles,(1-auxiliar.p.30)*60,col='orange', lwd=2) lines(horas.utiles,(1-auxiliar.p.20)*60,col='green', lwd=2) lines(horas.utiles,(1-auxiliar.p.10)*60,col='yellow', lwd=2) legend(7,60,legend=c('umbral = 40 km/h','umbral = 30 km/h','umbral = 20 km/h', 'umbral = 10 km/h'), col=c('red','orange','green','yellow'), lty=c(1,1)) # 4. Por último, se crea el vector con predicciones instantáneas hora.estudio <- as.data.frame(horas.utiles) predicciones.instantaneas <- (na.omit(prediccion.prueba$fit))^2 archivo.predicciones.instantaneas <- cbind(hora.estudio, predicciones.instantaneas)