srouchier / buildingenergygeeks

Bringing building energy simulation and statistical learning together
http://buildingenergygeeks.org
36 stars 15 forks source link

ARMAX Model in R using Rstan #4

Open bui090 opened 3 weeks ago

bui090 commented 3 weeks ago

Hi, I'm trying to build an ARMAX model as given in the book. Refer Chapter 7. I tried the tutorial as given: installed relevant packages, even installed Rtools as I was getting 'make' command error messages.

The code is running upto a certain extent. image

However, the issue I am facing is with the plotting of the graph using ggplot. The solver is not able to produce the graph and neither is it showing some error.

Plot

plotStart <- as.Date("2016-05-23", "%Y-%m-%d") # nolint # nolint plotEnd <- testEnd ggplot(data = df.test %>% filter(datetime > plotStart, datetime < plotEnd)) + geom_line(mapping = aes(x=ymd_hms(datetime), y=m0)) + geom_line(mapping = aes(x=ymd_hms(datetime), y=pred_med), color="red") + geom_ribbon(mapping = aes(x=ymd_hms(datetime), ymin=pred_low, ymax=pred_up), fill="red", alpha=0.1) ggsave("plot.png",width=5, [height=5)

I will appreciate any help if possible. Thanks

Find below the code I am using: library("readxl") library("languageserver") library("tidyverse") library("lubridate") library("dplyr") library("rstan") pkgbuild::has_build_tools(debug = TRUE)

df <- read.csv("downloads/building_1298.csv") head(df) df <- df %>% # nolint: trailing_whitespace_linter, trailing_whitespace_linter. mutate(Date = as_date(datetime), hour = hour(datetime), wday = wday(Date), week.end = wday==1 | wday==7) %>% fill(air_temperature, wind_speed)

trainingStart <- as.Date("2016-05-01", "%Y-%m-%d") trainingEnd <- as.Date("2016-05-31", "%Y-%m-%d") testEnd <- as.Date("2016-06-07", "%Y-%m-%d") df.train <- df %>% filter(datetime > trainingStart, datetime < trainingEnd) df.test <- df %>% filter(datetime > trainingStart, datetime < testEnd)

model = " data { int N_train; // data points (training) int N; // data points (all) int K; // number of predictors int P; // AR order

matrix[N, K] x; // predictor matrix (training and test) vector[N_train] y; // output data }

parameters { real alpha; // mean coeff row_vector[P] beta; // AR parameters row_vector[P+1] theta[K]; // coefficients of explanatory variables real sigma; // error term }

model { for (n in P+1:N_train) { real ex = 0; for (j in 1:K) ex += theta[j] x[n-P:n, j]; y[n] ~ normal(alpha + beta y[n-P:n-1] + ex, sigma); } }

generated quantities { vector[N] y_new;

for (n in 1:P) y_new[n] = y[n];

for (n in P+1:N_train) { real ex = 0; for (j in 1:K) ex += theta[j] x[n-P:n, j]; y_new[n] = normal_rng(alpha + beta y[n-P:n-1] + ex, sigma); }

for (n in N_train+1:N) { real ex = 0; for (j in 1:K) ex += theta[j] x[n-P:n, j]; y_new[n] = normal_rng(alpha + beta y_new[n-P:n-1] + ex, sigma); } } " model.data <- list( N_train = nrow(df.train), # number of training data points N = nrow(df.test), # total number of data points K = 1, # number of predictors P = 1, # AR order y = df.train$m0, # dependent variable (training data) x = df.test %>% select(air_temperature) # explanatory variables (all data) )

fit <- stan( model_code = model, # Stan program data = model.data, # named list of data chains = 4, # number of Markov chains warmup = 1000, # number of warmup iterations per chain iter = 2000, # total number of iterations per chain cores = 2, # number of cores (could use one per chain) refresh = 1, # progress not shown )

print(fit, pars=c("alpha", "beta", "theta", "sigma"))

la <- rstan::extract(fit, permuted = TRUE)

Quantiles from the posterior predictive distribution

y.quantiles <- apply(la$y_new, 2, quantile, probs=c(0.025, 0.5, 0.975)) df.test <- df.test %>% mutate(pred_low = y.quantiles[1,], pred_med = y.quantiles[2,], pred_up = y.quantiles[3,] )

Plot

plotStart <- as.Date("2016-05-23", "%Y-%m-%d") # nolint # nolint plotEnd <- testEnd ggplot(data = df.test %>% filter(datetime > plotStart, datetime < plotEnd)) + geom_line(mapping = aes(x=ymd_hms(datetime), y=m0)) + geom_line(mapping = aes(x=ymd_hms(datetime), y=pred_med), color="red") + geom_ribbon(mapping = aes(x=ymd_hms(datetime), ymin=pred_low, ymax=pred_up), fill="red", alpha=0.1) ggsave("plot.png",width=5, height=5)

srouchier commented 1 week ago

Hello I just ran the code from the ARMAX page without any issues. https://buildingenergygeeks.org/armax.html Could you please give more information on the errors you are having ?

bui090 commented 3 days ago

Hello I just ran the code from the ARMAX page without any issues. https://buildingenergygeeks.org/armax.html Could you please give more information on the errors you are having ?

I have shared the code I am using above. The issues is that despite making use of ggplot for plotting the graphs. I just get a blank graph. I have ensured I have the proper libraries installed, made sure no semantics error is there in the code. There is no explicit error message. I do not know is its an error with R or visual studio? Appreciate any assistance. Thanks