CambridgeTemperatureModel

Time series and other models for Cambridge UK temperature forecasts in R

If you like CambridgeTemperatureModel, give it a star, or fork it and contribute!

Summary of forecast error for mean, naive, simple exponential smoothing and Holt smoothing methods for horizons up to 24 hours (48 * 1/2 hour forecasts):

Further details in Daily forecast baselines section below.

NOTE: Data set has been substantially updated and cleaned since graph creation.

Installation/Usage

Requires R version 3.2.0 and higher.

To install the required libraries in an R session:

install.packages("caret", dependencies = c("Depends", "Suggests"))
install.packages("data.table")
install.packages("stationaRy")
install.packages("lubridate")
install.packages("forecast")
install.packages("prophet")
install.packages("suncalc")
install.packages("ggplot2")
install.packages("tseries")

Clone repository:

git clone https://github.com/makeyourownmaker/CambridgeTemperatureModel
cd CambridgeTemperatureModel

The R files can be ran in sequence or the R session image can be loaded.

Run files in sequence in an R session:

setwd("CambridgeTemperatureModel")
source("1-load.R", echo = TRUE)
source("2-clean.R", echo = TRUE)

Or load R session image:

setwd("CambridgeTemperatureModel")
load("data/CambridgeTemperatureModel.RData")

Details

The Digital Technology Group in the Cambridge University Computer Laboratory maintain a weather station.

I live close to this weather station. When I started looking at this data the UK Met Office were updating forecasts every 2 hours. I thought I could produce a more frequent nowcast (one step ahead forecast) using time series or statistical learning methods. Day long forecasts are of secondary interest. Temperature and rainfall are the primary variables of interest. Unfortunately, the rain sensor has issues.

I have no affiliation with Cambridge University, the Computer Laboratory or the Digital Technology Group.

Variables

The weather measurements include the following variables.

Variables	Units
Temperature	Celsius (°C) * 10
Dew Point	Celsius (°C) * 10
Humidity	Percent
Pressure	mBar
Wind Speed Mean	Knots * 10
Wind Bearing Mean	Degrees
Timestamp	Data Hours:Minutes:Seconds

Dew point is the temperature at which air, at a level of constant pressure, can no longer hold all the water it contains. Dew point is defined here and in more detail here.

There are known issues with the sunlight and rain sensors. These measurements are not included for now.

Measurements are recorded every 30 minutes.

Cleaning

The data included start on 2008-08-01 when the weather station was moved to it's current location.

Unfortunately, the data is quite noisy and usually have a couple of hundred missing observations every year. The following cleaning steps are performed:

• The Digital Technology Group list inaccuracies in the weather data which were excluded.
• Unrealistically high and low values were removed.
• Long runs of consecutively equal values were eliminated.
• Sudden large increasing/decreasing observations were omitted.
• Cook's distance was used to remove influential observations.
• Cambridge Airport weather measurements from ISD were used to find additional outliers in the Computer Lab measurements.
  • the above exclusions were also applied to the ISD data
• Limited (12 hours max) linear interpolation was used to fill missing observations.
• Where available Cambridge Airport ISD data was used to fill missing/deleted values.
  • no pressure data for Cambridge Airport from ISD
• Historical averages were used to fill missing pressure, wind speed and bearing values.
  • multiple imputation for these variables gave poor results
• Multiple imputation was used to replace missing temperature, humidity and dew point values.
  • multiple imputation for these variables gave similar or better results to historical averages

The most recent cleaned data have no missing values. Data older than 2021/04/26 have had less cleaning. Outlier exclusion has been fairly conservative. Some problems may remain in the data, such as short and/or long term sensor drift or periods of anomolously high variance.

The following figure shows an older cleaned temperature time series.

A visual inspection indicates a lack of trend.

The ADF and KPSS tests in the exloratory data analysis file (described in the Files subsection below) implies the stationarity of this time series.

Cambridge Airport ISD data

Cambridge Airport weather measurements from ISD were used to find outliers in the Computer Laboratory measurements and to replace missing values. The stationaRy package was used to download the ISD data. Unfortunately there are no pressure measurements in the Airport observations. The ISD data is somewhat cleaner than the Computer lab data. Data cleaning and limited interpolation were applied to the Cambridge Airport data before being used to replace NAs in the Computer lab data.

One step ahead baselines

The following table shows accuracy metrics for baseline nowcast methods:

Method	RMSE	MAE	MAPE
Mean temperature	64.46	52.63	249.91
Persistent temperature	6.26	4.13	9.49
Simple exponential smoothing	6.05	4.03	9.81
Holt exponential smoothing	5.62	3.94	10.25

These results are from older partially cleaned observations.

These metrics are calculated in the baselines file briefly described in the Files subsection. Numbers in bold indicate the lowest value for each metric.

The three accuracy metrics:

• RMSE - Root Mean Squared Error
• MAE - Mean Absolute Error
• MAPE - Mean Absolute Percent Error

The four baseline methods:

1. The mean temperature method simply uses the mean temperature across the entire data set as the nowcast.
2. The persistent temperature method is a popular benchmark in the meteorology literature. It uses the previous temperature value for the nowcast. The forecast package documentation refers to it as the naive method.
3. Simple exponential smoothing (ses) uses "weighted averages, where the weights decrease exponentially as observations come from further in the past". Generally speaking, this method is surprisingly accurate given its low computational complexity.
4. Holt extended simple exponential smoothing to include data with a trend.

Daily forecast baselines

The following graph shows RMSE values for baseline daily forecast methods:

The ses and naive methods give almost identical results.

Two different Holt-Winters exponential smoothing implementations failed! Sadly the double seasonal Holt-Winters exponential smoothing implementation in the forecast package is not suitable when data contain zeros or negative numbers. Vanilla ARIMA models are not suitable for this temperature data due to multi-seasonality which is explained next.

Seasonality

In general, time series can be decomposed into seasonal and trend components.

The Cambridge temperature data contains two seasonal components:

1. Daily
2. Yearly

The next two figures show the daily and yearly components found using the prophet package. This code is briefly described in the Files subsection.

1. Daily seasonal trend component

2. Yearly seasonal trend component

The daily and yearly components show smooth cyclic change as expected. The vertical axis shows percent change in temperature.

Prophet package models

Prophet models are robust to missing data, shifts in the trend and typically handle outliers well. Yearly, weekly, and daily seasonality, plus holiday effects can be accomodated. Seasonal components are represented using Fourier terms. Prophet models work best with time series that have strong seasonal effects and several seasons of historical data. Stan is used for fitting models.

Two prophet models were built:

1. Logistic growth with daily and annual components with automatic changepoint detection
2. Logistic growth with daily and annual components with 50 changepoints specified

In both cases a floor of -150 and a cap of 400 were used for logistic growth.

A changepoint is a specific timepoint where the statistical properties differ before and after the timepoint. The prophet package detects 25 changepoints automatically.

Additive seasonality is assumed for both models.

The accuracy results for one step ahead forecasts:

Method	RMSE	MAE	MAPE
Logistic growth, automatic changepoints	28.82	25.88	50.25
Logistic growth, 50 changepoints	28.66	25.80	50.13

These results are from older partially cleaned observations.

Using more changepoints showed little to no improvement.

These results are substantially higher than most of the baseline one step ahead forecasts. It's possible that using more data would improve the yearly seasonal component and in turn improve the nowcasts. The prophet models may perform better for daily forecasts. Unfortunately, daily forecast cross-validation will be quite time-consuming to run.

Forecast package model

The forecast package supports multi-seasonal models using the TBATS (Trigonometric Exponential Smoothing) method.

This function uses a trigonometric representation of seasonality, instead of conventional seasonal indices. It also automatically performs Box-Cox transformation of the time series, as required. It can be very slow to estimate, especially with multiple seasonal time series. The tbats() function does not support including additional regressors.

Unfortunately, cross-validation fails. See the source code described in the Files subsection for details and this unanswered stackoverflow question.

FWIW here are the training set accuracy metrics for one step ahead forecasts:

Method	RMSE	MAE	MAPE
TBATS	5.7	3.8	Inf

These results are from older partially cleaned observations.

These results are not comparable with the baseline methods which are calculated on a separate test data set.

The infinite MAPE value comes from the forecast package mape() function implementation which permits division by zero. Other implementations add one to the denominator to avoid this behavior.

Files

These files demonstrate how to build models for the Cambridge UK temperature data:

• 1-load.R
  • Download data, set variable types and adds some date and time related fields
    • Both computer lab and NOAA ISD Cambridge airport data
• 2-clean.R
  • Remove known inaccuracies and other unrealistic measurements
  • Use NOAA ISD Cambridge airport data, historical averages and multiple imputation to replace missing values
• 3.01-eda.R
  • Temperature time series plot and stationarity tests
  • This script will create a directory called figures if it doesn't already exist
  • Exploratory data analysis is also covered in my Cambridge University Computer Laboratory Weather Station R Shiny Web App repository
• Some feature engineering will be required
  • Transformations like the Box-Cox
  • Dummy seasonal variables for certain models
  • Possibly deseasonalisation
• 4.01-baselines.R
  • Build baseline models and calculate nowcast and daily accuracy using the forecast package.
    • This script will create a directory called figures if it doesn't already exist
• 4.02-prophet.R
  • Build multi-seasonal model using the prophet package.
    • This script will create a directory called figures if it doesn't already exist
• 4.03-forecast.R
  • Build multi-seasonal TBATS model using the forecast package.
    • Cross-validation fails - see source code for details

Roadmap

• Add standard deviations to MSE, MAE and MAPE values
• Further develop data cleaning
  • Cook's distance based outlier detection
    • Would benefit from a seasonal component (loess smoothed 99th percentile) instead of single static value
    • Consider re-running for individual features and/or sub-groups of features
  • Cook's distance alternatives for finding and replacing outliers
    • tsrobprep Robust Preprocessing of Time Series Data
      • Particularly the robust time series seasonal decomposition function robust_decompose
    • ctbi A Procedure to Clean, Decompose and Aggregate Timeseries
    • rucm Implementation of Unobserved Components Model
      • Supports cycles (eg annual seasonality) and shorter (eg daily) seasonalities plus regressors
• Consider alternatives to stationaRy
  • stationaRy supports years upto 2020 so requires monkey-patching the get_years_available_for_station function
  • Investigate missing airport pressure data from stationaRy
  • worldmet Import Surface Meteorological Data from NOAA Integrated Surface Database (ISD)
  • ropensci/riem queries global ASOS data from IEM
    • ASOS data from IEM has 30 min updates for Cambridge airport
  • ropensci/rnoaa Client for many 'NOAA' data sources
  • ropensci/nasapower API Client for NASA POWER Global Meteorology, Surface Solar Energy and Climatology
    • Satellite observations and calculations
    • Updated daily but main features (temperature, dew.point, pressure, humidity) not available for previous 3 days
    • Additional features (longwave radiation, shortwave radiation etc) have longer latencies
    • Potentially useful for replacing long sequences of multiple missing features
• Examine Global Forecast System (GFS) weather model
  • Runs four times a day, produces forecasts up to 16 days in advance
  • Data is available for free in the public domain
  • Model serves as the basis for the forecasts of numerous services
  • Potentially use for
    • additional exogeneous variables and/or future covariates
    • model averaging
    • interesting benchmark if nothing else
  • Unsure about GFS forecast accuracy - could be adding more variance
    • Probably not for filling missing values
  • Don't want to directly deal with GRIB or netcdf files
    • meteoForecast access real time and historic Numerical Weather Predictions as a time series for a location
    • rNOMADS An interface to the NOAA Operational Model Archive and Distribution System (NOMADS)
      • List of NOMADS global models
• Feature calculations
  • humidity Calculate Water Vapor Measures from Temperature and Dew Point
  • Ground heat flux
  • Calculate longwave and shortwave radiation approximations
    • Evaluation of Clear-Sky Incoming Radiation Estimating Equations Typically Used in Remote Sensing Evapotranspiration Algorithms
      • Clear-sky only :-(
• Investigate most extreme changepoints in observations
  • Such as short term sensor drift or high variance periods
  • Using strucchange or a similar package
• Add more time series models
  • I have some GAM forecasts which are nearing completion
  • TSA supports multiple seasonalities and exogenous variables with the arimax() function
• Improve documentation
  • Describe cross-validation

Contributing

Pull requests are welcome. For major changes, please open an issue first to discuss what you would like to change.

Alternatives

• UK Met Office
• Cambridge University Computer Laboratory Weather Station R Shiny Web App
• Forecasting surface temperature based on latitude, longitude, day of year and hour of day
• Time series and other models for Cambridge UK temperature forecasts in python

License

GPL-2