# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     name: python3
# ---

# + [markdown] id="view-in-github" colab_type="text"
# <a href="https://colab.research.google.com/github/desireechen/Deep-Learning/blob/master/Time_series%2C_sequences_and_predictions.ipynb" target="_parent"><img src="https://colab.research.google.com/assets/colab-badge.svg" alt="Open In Colab"/></a>

# + id="aCMZXB22WjWm" colab_type="code" colab={}
import numpy as np
import matplotlib.pyplot as plt
import tensorflow as tf
from tensorflow import keras


# + id="PEQRASQ43c9E" colab_type="code" colab={}
# Get sunspot data.
# !wget --no-check-certificate \
#     https://storage.googleapis.com/laurencemoroney-blog.appspot.com/Sunspots.csv \
#     -O /tmp/sunspots.csv

# + id="TTOrgtv7Cao5" colab_type="code" colab={}
# Get temperature data.
# !wget --no-check-certificate \
#     https://raw.githubusercontent.com/jbrownlee/Datasets/master/daily-min-temperatures.csv \
#     -O /tmp/daily-min-temperatures.csv

# + id="3nc2ss9BWk2B" colab_type="code" colab={}
# OPTION 1 for plot series function.
def plot_series(time, series, format="-", start=0, end=None, label=None): # Sometimes do not have label argument.
    plt.plot(time[start:end], series[start:end], format, label=label) # Sometimes do not have label.
    plt.xlabel("Time")
    plt.ylabel("Value")
    if label:
        plt.legend(fontsize=14)
    plt.grid(True) # Can use either True or False.


# + id="P_zxv13SjFih" colab_type="code" colab={}
# OPTION 2 for plot series function.
def plot_series(time, series, format="-", start=0, end=None):
    plt.plot(time[start:end], series[start:end], format)
    plt.xlabel("Time")
    plt.ylabel("Value")
    plt.grid(True) # Can use either True or False.


# + id="km3XAnWg5HEs" colab_type="code" colab={}
# Process sunspot data.
import csv
time_step = []
sunspots = []

with open('/tmp/sunspots.csv') as csvfile:
  reader = csv.reader(csvfile, delimiter=',')
  next(reader)
  for row in reader:
    time_step.append(int(row[0])) #1st column contains the index number.
    sunspots.append(float(row[2])) #3rd column contains the sunspot.
# Make numpy array.
time = np.array(time_step)
series = np.array(sunspots)
# Plot.
plt.figure(figsize=(10, 6))
plot_series(time, series) # To see both train and test. Have not split yet.
plt.show()

# + id="G3khYEZNC2gw" colab_type="code" colab={}
# Process temperature data.
import csv
time_step = []
temps = []

with open('/tmp/daily-min-temperatures.csv') as csvfile:
  reader = csv.reader(csvfile, delimiter=',')
  next(reader)
  step=0
  for row in reader:
    time_step.append(step)
    temps.append(float(row[1]))
    step = step + 1
# Make numpy array.
time = np.array(time_step)
series = np.array(temps)
# Plot.
plt.figure(figsize=(10, 6))
plot_series(time, series) # To see both train and test. Have not split yet.
plt.show()


# + id="mB9qtcBeWwP5" colab_type="code" colab={}
def trend(time, slope=0):
    return slope * time


# + id="ro4vh3SaXJlG" colab_type="code" colab={}
# OPTION 1 for seasonal pattern function.
def seasonal_pattern(season_time):
    """Just an arbitrary pattern, you can change it if you wish"""
    return np.where(season_time < 0.4,
                    np.cos(season_time * 2 * np.pi),
                    1 / np.exp(3 * season_time))


# + id="wjJ19KwFiDdk" colab_type="code" colab={}
# OPTION 2 for seasonal pattern function.
def seasonal_pattern(season_time):
    """Just an arbitrary pattern, you can change it if you wish"""
    return np.where(season_time < 0.1,
                    np.cos(season_time * 7 * np.pi),
                    1 / np.exp(5 * season_time))


# + id="a8zOV9TUy9pD" colab_type="code" colab={}
# OPTION 3 for seasonal pattern function.
def seasonal_pattern(season_time):
    """Just an arbitrary pattern, you can change it if you wish"""
    return np.where(season_time < 0.1,
                    np.cos(season_time * 6 * np.pi),
                    2 / np.exp(9 * season_time))


# + id="2oM03z5SXKF-" colab_type="code" colab={}
def seasonality(time, period, amplitude=1, phase=0):
    """Repeats the same pattern at each period"""
    season_time = ((time + phase) % period) / period
    return amplitude * seasonal_pattern(season_time)


# + id="wN7gEWXUXd4g" colab_type="code" colab={}
def noise(time, noise_level=1, seed=None):
    rnd = np.random.RandomState(seed)
    return rnd.randn(len(time)) * noise_level


# + id="DFMFZxBFlvl9" colab_type="code" colab={}
# For no Conv1D. This is to turn a series into a dataset which we can train on.
def windowed_dataset(series, window_size, batch_size, shuffle_buffer_size):
  dataset = tf.data.Dataset.from_tensor_slices(series)
  dataset = dataset.window(window_size + 1, shift=1, drop_remainder=True)
  dataset = dataset.flat_map(lambda window: window.batch(window_size + 1))
  dataset = dataset.shuffle(shuffle_buffer_size).map(lambda window: (window[:-1], window[-1]))
  dataset = dataset.batch(batch_size).prefetch(1)
  return dataset


# + id="b5SKzG8S69F_" colab_type="code" colab={}
# For Conv1D. This is to turn a series into a dataset which we can train on.
def windowed_dataset(series, window_size, batch_size, shuffle_buffer_size):
  series = tf.expand_dims(series, axis=-1) 
  dataset = tf.data.Dataset.from_tensor_slices(series)
  dataset = dataset.window(window_size + 1, shift=1, drop_remainder=True)
  dataset = dataset.flat_map(lambda w: w.batch(window_size + 1))
  dataset = dataset.shuffle(shuffle_buffer_size)
  dataset = dataset.map(lambda w: (w[:-1], w[1:]))
  return dataset.batch(batch_size).prefetch(1)


# + id="nx71ecypwH2r" colab_type="code" colab={}
def model_forecast(model, series, window_size):
    ds = tf.data.Dataset.from_tensor_slices(series)
    ds = ds.window(window_size, shift=1, drop_remainder=True)
    ds = ds.flat_map(lambda w: w.batch(window_size))
    ds = ds.batch(32).prefetch(1)
    forecast = model.predict(ds)
    return forecast


# + id="ZUCshmjveQjw" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 388} outputId="009b18fd-5cbe-4be4-b620-b94d48c31456"
time = np.arange(4 * 365 + 1, dtype="float32") # Instead of 4, it can be 10.
series = trend(time, 0.1)  
baseline = 10
amplitude = 40 # This can be 20 or 40.
slope = 0.01 # This can be 0.01, 0.05, 0.09 or 0.005.
noise_level = 2 # This can be 2, 3 or 5.
# Create the series.
series = baseline + trend(time, slope) + seasonality(time, period=365, amplitude=amplitude)
# Update with noise.
series += noise(time, noise_level, seed=42) # Can try other seed such as 51.
# Plot.
plt.figure(figsize=(10, 6))
plot_series(time, series) # To see both train and test. Have not split yet.
plt.show()

# + id="Vt90mG-vX8Oy" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 759} outputId="3d1e7759-6441-4b2e-cf5e-3e20b1a360a2"
split_time = 1100 # This number can be 1000, 1100, 2500 or 3000.
time_train = time[:split_time]
x_train = series[:split_time]
time_valid = time[split_time:]
x_valid = series[split_time:]

plt.figure(figsize=(10, 6))
plot_series(time_train, x_train)
plt.show()
plt.figure(figsize=(10, 6))
plot_series(time_valid, x_valid)
plt.show()

# + id="2iY3gyERmZdg" colab_type="code" colab={}
# Always clear session before training another model.
tf.keras.backend.clear_session() # Clears any internal variables.
tf.random.set_seed(51)
np.random.seed(51)

# + id="KUVtzhLImQYJ" colab_type="code" colab={}
# Windowed dataset.
window_size = 20 # This can be 20, 30, 60 or 64.
batch_size = 32 # This can be 32, 100, 128 or 256.
shuffle_buffer_size = 1000

dataset = windowed_dataset(x_train, window_size, batch_size, shuffle_buffer_size)
# print(dataset)
# print(x_train.shape)

# + id="pDOwATxAsl2S" colab_type="code" colab={}
# Single layer.
single_layer = tf.keras.layers.Dense(1, input_shape=[window_size])
model = tf.keras.models.Sequential([single_layer])
# model.summary()
model.compile(optimizer=tf.keras.optimizers.SGD(lr=1e-6, momentum=0.9), loss='mse')

print('Layer weights {}'.format(single_layer.get_weights())) # 1st array has 20 weights for the 20 x values.
                                                            # 2nd array is the Bias.

print(series[1:21])
model.predict(series[1:21][np.newaxis])

# + id="HzFX6Kaxuy80" colab_type="code" colab={}
# DNN with Learning Rate Scheduler.
model = tf.keras.models.Sequential([tf.keras.layers.Dense(10, input_shape=[window_size], activation='relu'), 
                                    # Number above can be 10, 20 or 100.
                                    tf.keras.layers.Dense(10, activation='relu'),
                                    tf.keras.layers.Dense(1)])
# model.summary()
# Learning Rate Scheduler will adjust the learning rate over each epoch.
lr_schedule = tf.keras.callbacks.LearningRateScheduler(
    lambda epoch: 1e-8 * 10**(epoch / 20)) # LR same as below and vice versa. Can be lr=1e-6, 1e-8, 8e-6, 5e-5, 1e-5 or 1e-7.
model.compile(optimizer=tf.keras.optimizers.SGD(lr=1e-8, momentum=0.9), loss='mse', metrics=['mae']) 
num_epochs = 100 # This number can be 10, 30, 50, 100, 150, 400 or 500.
history = model.fit(dataset, epochs=num_epochs, callbacks=[lr_schedule], verbose=0)
# verbose 0, 1, 2 shows the least, most, in between information.

# + id="D6CEL7qmsnsi" colab_type="code" colab={}
# DNN without Learning Rate Scheduler.
model = tf.keras.models.Sequential([tf.keras.layers.Dense(10, input_shape=[window_size], activation='relu'), 
                                    # Number above can be 10, 20 or 100.
                                    tf.keras.layers.Dense(10, activation='relu'),
                                    tf.keras.layers.Dense(1)])
# model.summary()
model.compile(optimizer=tf.keras.optimizers.SGD(lr=1e-6, momentum=0.9), loss='mse', metrics=['mae']) # Can be lr=1e-6, 1e-8, 8e-6, 5e-5, 1e-5 or 1e-7.
num_epochs = 100 # This number can be 10, 30, 50, 100, 150, 400 or 500.
model.fit(dataset, epochs=num_epochs,verbose=0)
# verbose 0, 1, 2 shows the least, most, in between information.

# + id="0K1ajFkvm-ob" colab_type="code" colab={}
# RNN with Learning Rate Scheduler.
model = tf.keras.models.Sequential([tf.keras.layers.Lambda(lambda x: tf.expand_dims(x, axis=-1), input_shape=[None]), 
                                    tf.keras.layers.SimpleRNN(40, return_sequences=True),
                                    tf.keras.layers.SimpleRNN(40),
                                    tf.keras.layers.Dense(1),
                                    tf.keras.layers.Lambda(lambda x: x*100.0)]) # This number can be 100.0 or 200 or 400, depends on range of numbers.
# Learning Rate Scheduler will adjust the learning rate over each epoch.
lr_schedule = tf.keras.callbacks.LearningRateScheduler(
    lambda epoch: 1e-8 * 10**(epoch / 20)) # LR same as below and vice versa. Can be lr=1e-6, 1e-8, 8e-6, 5e-5, 1e-5 or 1e-7.
model.compile(optimizer=tf.keras.optimizers.SGD(lr=1e-8, momentum=0.9), loss=tf.keras.losses.Huber(), metrics=['mae']) 
num_epochs = 100 # This number can be 10, 30, 50, 100, 150, 400 or 500.
history = model.fit(dataset, epochs=num_epochs,callbacks=[lr_schedule]) # Did not include verbose argument, so default.
# verbose 0, 1, 2 shows the least, most, in between information.

# + id="ol7wA3JJs-vb" colab_type="code" colab={}
# RNN without Learning Rate Scheduler.
model = tf.keras.models.Sequential([tf.keras.layers.Lambda(lambda x: tf.expand_dims(x, axis=-1), input_shape=[None]), 
                                    tf.keras.layers.SimpleRNN(40, return_sequences=True),
                                    tf.keras.layers.SimpleRNN(40),
                                    tf.keras.layers.Dense(1),
                                    tf.keras.layers.Lambda(lambda x: x*100.0)]) # This number can be 100.0 or 200 or 400, depends on range of numbers.
model.compile(optimizer=tf.keras.optimizers.SGD(lr=5e-5, momentum=0.9), loss=tf.keras.losses.Huber(), metrics=['mae'])
                                                  # Can be lr=1e-6, 1e-8, 8e-6, 5e-5, 1e-5 or 1e-7.
num_epochs = 100 # This number can be 10, 30, 50, 100, 150, 400 or 500.
history = model.fit(dataset, epochs=num_epochs) # Did not include verbose argument.
# verbose 0, 1, 2 shows the least, most, in between information.

# + id="7owCwFH0wrIH" colab_type="code" colab={}
# LSTM with Learning Rate Scheduler.
model = tf.keras.models.Sequential([tf.keras.layers.Lambda(lambda x: tf.expand_dims(x, axis=-1), input_shape=[None]),
                                    tf.keras.layers.Bidirectional(tf.keras.layers.LSTM(32, return_sequences=True)),
                                    tf.keras.layers.Bidirectional(tf.keras.layers.LSTM(32)),
                                    # LSTM number can be 32, 60 or 64.
                                    tf.keras.layers.Dense(1),
                                    tf.keras.layers.Lambda(lambda x: x*100.0)]) # This number can be 100.0 or 200 or 400, depends on range of numbers.
# Learning Rate Scheduler will adjust the learning rate over each epoch.
lr_schedule = tf.keras.callbacks.LearningRateScheduler(
    lambda epoch: 1e-8 * 10**(epoch / 20)) # LR same as below and vice versa. Can be lr=1e-6, 1e-8, 8e-6, 5e-5, 1e-5 or 1e-7.
model.compile(optimizer=tf.keras.optimizers.SGD(lr=1e-8, momentum=0.9), loss=tf.keras.losses.Huber(), metrics=['mae']) 
num_epochs = 100 # This number can be 10, 30, 50, 100, 150, 400 or 500.
history = model.fit(dataset, epochs=num_epochs, callbacks=[lr_schedule]) # Did not include verbose argument.
# verbose 0, 1, 2 shows the least, most, in between information.

# + id="cqJuwW1H0rGf" colab_type="code" colab={}
# LSTM without Learning Rate Scheduler. Can have multi-layer LSTMs.
model = tf.keras.models.Sequential([tf.keras.layers.Lambda(lambda x: tf.expand_dims(x, axis=-1), input_shape=[None]),
                                    tf.keras.layers.Bidirectional(tf.keras.layers.LSTM(32, return_sequences=True)),
                                    # tf.keras.layers.Bidirectional(tf.keras.layers.LSTM(32, return_sequences=True)),
                                    tf.keras.layers.Bidirectional(tf.keras.layers.LSTM(32)),
                                    # LSTM number can be 32, 60 or 64.
                                    tf.keras.layers.Dense(1),
                                    tf.keras.layers.Lambda(lambda x: x*100.0)]) # This number can be 100.0 or 200 or 400, depends on range of numbers.
model.compile(optimizer=tf.keras.optimizers.SGD(lr=1e-5, momentum=0.9), loss='mse', metrics=['mae']) 
                                          # Sometimes do not have metrics argument. Then, it is not shown in output.
                                          # Can be lr=1e-6, 1e-8, 8e-6, 5e-5, 1e-5 or 1e-7.
num_epochs = 100 # This number can be 10, 30, 50, 100, 150, 400 or 500.
history = model.fit(dataset, epochs=num_epochs) # Did not include verbose argument.
# verbose 0, 1, 2 shows the least, most, in between information.

# + id="BAoZ7UCIyigg" colab_type="code" colab={}
# Conv1D and LSTM with Learning Rate Scheduler.
model = tf.keras.models.Sequential([tf.keras.layers.Conv1D(filters=32, kernel_size=5, strides=1,
                                                           padding='causal', activation='relu', input_shape=[None, 1]),
                                    # filters can be 32 or 60. kernel size 3 or 5. 
                                    tf.keras.layers.Bidirectional(tf.keras.layers.LSTM(32, return_sequences=True)), 
                                    # LSTM number can be 32, 60 or 64.
                                    tf.keras.layers.Bidirectional(tf.keras.layers.LSTM(32, return_sequences=True)),
                                    # tf.keras.layers.Dense(30, activation='relu'),
                                    # tf.keras.layers.Dense(10, activation='relu'),
                                    tf.keras.layers.Dense(1),
                                    tf.keras.layers.Lambda(lambda x: x*100.0)]) # This number can be 100.0 or 200 or 400, depends on range of numbers.
# Learning Rate Scheduler will adjust the learning rate over each epoch.
lr_schedule = tf.keras.callbacks.LearningRateScheduler(
    lambda epoch: 1e-8 * 10**(epoch / 20)) # LR same as below and vice versa. Can be lr=1e-6, 1e-8, 8e-6, 5e-5, 1e-5 or 1e-7.
model.compile(optimizer=tf.keras.optimizers.SGD(lr=1e-8, momentum=0.9), loss=tf.keras.losses.Huber(), metrics=['mae']) 
num_epochs = 100 # This number can be 10, 30, 50, 100, 150, 400 or 500.
history = model.fit(dataset, epochs=num_epochs, callbacks=[lr_schedule]) # Did not include verbose argument.
# verbose 0, 1, 2 shows the least, most, in between information.

# + id="CaMcP44Q0Z4s" colab_type="code" colab={}
# Conv1D and LSTM without Learning Rate Scheduler.
model = tf.keras.models.Sequential([tf.keras.layers.Conv1D(filters=32, kernel_size=5, strides=1, 
                                                           padding='causal', activation='relu', input_shape=[None, 1]),
                                    # filters can be 32 or 60. kernel size 3 or 5. 
                                    tf.keras.layers.LSTM(32, return_sequences=True), # LSTM number can be 32, 60 or 64.
                                    tf.keras.layers.LSTM(32, return_sequences=True), # LSTM number can be 32, 60 or 64.
                                    # tf.keras.layers.Dense(30, activation='relu'),
                                    # tf.keras.layers.Dense(10, activation='relu'),
                                    tf.keras.layers.Dense(1),
                                    tf.keras.layers.Lambda(lambda x: x*100.0)]) # This number can be 100.0 or 200 or 400, depends on range of numbers.
model.compile(optimizer=tf.keras.optimizers.SGD(lr=1e-8, momentum=0.9), loss=tf.keras.losses.Huber(), metrics=['mae'])
                                                # Can be lr=1e-6, 1e-8, 8e-6, 5e-5, 1e-5 or 1e-7.
num_epochs = 100 # This number can be 10, 30, 50, 100, 150, 400 or 500.
history = model.fit(dataset, epochs=num_epochs) # Did not include verbose argument.
# verbose 0, 1, 2 shows the least, most, in between information.

# + id="UZ1q8Or7xCol" colab_type="code" colab={}
# Extra code. Plot loss against learning rates. From plot, choose learning rate with minimum loss.
## EITHER THIS.
lrs = 1e-8 * (10 ** (np.arange(100) / 20))
## OR THIS.
lrs = history.history['lr']
losses = history.history['loss']
plt.semilogx(lrs, losses)
plt.axis([1e-8, 1e-3, 0, 300]) # Can be [1e-8, 1e-4, 0, 30] or [1e-8, 1e-4, 0, 60]
                              # 1st 2 is the range of x-axis, next 2 is the range of the y-axis.

# + id="mN-fSaMYxzgM" colab_type="code" colab={}
# Extra code. Plot loss against epoch.
loss = history.history['loss']
epochs = range(len(loss))
plt.plot(epochs, loss, 'b', label='Training Loss')
plt.show()

# Remove the 1st 10 epochs.
loss = history.history['loss']
epochs = range(10, len(loss))
plot_loss = loss[10:]
print(plot_loss)
plt.plot(epochs, plot_loss, 'b', label='Training Loss')
plt.show()

# + id="509RnqRboeH4" colab_type="code" colab={}
# Forecasting. EITHER THIS.
forecast = []
results = [] 

for time in range(len(series) - window_size):
  forecast.append(model.predict(series[time:time + window_size][np.newaxis]))

forecast = forecast[split_time-window_size:] # Get the forecast after the split.
results = np.array(forecast)[:, 0, 0]

plt.figure(figsize=(10, 6))
plot_series(time_valid, x_valid)
plot_series(time_valid, results)

# Print metrics MAE and MSE.
print(keras.metrics.mean_absolute_error(x_valid, results).numpy())
print(tf.keras.metrics.mean_absolute_error(x_valid, results).numpy()) # If did not use from tensorflow import keras
print(keras.metrics.mean_squared_error(x_valid, results).numpy())
print(tf.keras.metrics.mean_squared_error(x_valid, results).numpy()) # If did not use from tensorflow import keras

# + id="cU1MCJnz1bX9" colab_type="code" colab={}
# Forecasting. OR THIS.
rnn_forecast = model_forecast(model, series[..., np.newaxis], window_size)
rnn_forecast = rnn_forecast[split_time - window_size:-1, -1, 0]
# print(rnn_forecast)

plt.figure(figsize=(10, 6))
plot_series(time_valid, x_valid)
plot_series(time_valid, rnn_forecast)

# Print metrics MAE and MSE.
print(keras.metrics.mean_absolute_error(x_valid, rnn_forecast).numpy())
print(tf.keras.metrics.mean_absolute_error(x_valid, rnn_forecast).numpy()) # If did not use from tensorflow import keras
print(keras.metrics.mean_squared_error(x_valid, rnn_forecast).numpy())
print(tf.keras.metrics.mean_squared_error(x_valid, rnn_forecast).numpy()) # If did not use from tensorflow import keras

# + id="jQKJg70vuIlG" colab_type="code" colab={}
import matplotlib.pyplot as plt # If this is already imported before, then do not need to import again. 
# import matplotlib.image as mpimg

mae=history.history['mae']
loss=history.history['loss']

epochs=range(len(loss)) # Get the number of epochs.

plt.plot(epochs, mae, 'r')
plt.plot(epochs, loss, 'b')
plt.title('MAE and Loss')
plt.xlabel("Epochs")
plt.ylabel("Accuracy")
plt.legend(["MAE", "Loss"])
plt.figure()

# Zoomed-in plot.
epochs_zoom = epochs[200:] # OR epochs_zoom = range(200,500)
mae_zoom = mae[200:]
loss_zoom = loss[200:]

plt.plot(epochs_zoom, mae_zoom, 'r')
plt.plot(epochs_zoom, loss_zoom, 'b')
plt.title('MAE and Loss')
plt.xlabel("Epochs")
plt.ylabel("Accuracy")
plt.legend(["MAE", "Loss"])
plt.figure()
# If see small spikes in the graphs, it may mean that the batch size is too small. Try using a larger batch size. 

# + id="ryqiZAidkInS" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 388} outputId="51a229a0-f6ca-45e1-f42c-24f28a55b0e8"
naive_forecast = series[split_time - 1:-1]
plt.figure(figsize=(10, 6))
plot_series(time_valid, x_valid)
plot_series(time_valid, naive_forecast)

# + id="8qVPFiMCkZLE" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 388} outputId="d9df8d0d-db65-47c7-b94e-8906601c121b"
plt.figure(figsize=(10, 6))
plot_series(time_valid, x_valid, start=0, end=150)
plot_series(time_valid, naive_forecast, start=1, end=151)

# + id="QHTgOn3Uk8Cl" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 52} outputId="3f8988db-ec39-4093-fdff-6bdaa127bcd9"
# Compute MSE and MAE between forecasts and predictions in the validation period. This will be the baseline. 
print(keras.metrics.mean_squared_error(x_valid, naive_forecast).numpy())
print(keras.metrics.mean_absolute_error(x_valid, naive_forecast).numpy())


# + id="hnx0rTDQfLmK" colab_type="code" colab={}
def moving_average_forecast(series, window_size):
  """Forecasts the mean of the last few values.
     If window_size=1, then this is equivalent to naive forecast"""
  forecast = []
  for time in range(len(series) - window_size):
    forecast.append(series[time:time + window_size].mean())
  return np.array(forecast)


# + id="SM6R7vQblO1V" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 388} outputId="6ec88354-06c2-47e2-8788-ea6d1f43fdbd"
moving_avg = moving_average_forecast(series, 30)[split_time - 30:]
plt.figure(figsize=(10, 6))
plot_series(time_valid, x_valid)
plot_series(time_valid, moving_avg)

# + id="KqGzedyMlVNX" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 52} outputId="2d354014-5522-4f04-fb54-201cd1910030"
# Compute MSE and MAE between Moving Average forecasts and predictions in the validation period. This Moving Average is window size 30. 
print(keras.metrics.mean_squared_error(x_valid, moving_avg).numpy())
print(keras.metrics.mean_absolute_error(x_valid, moving_avg).numpy())

# + id="6nVf8owNlt5W" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 388} outputId="0d044e43-36c5-4a94-91b4-7c2f78b9bfa5"
diff_series = (series[365:] - series[:-365])
diff_time = time[365:]
plt.figure(figsize=(10, 6))
plot_series(diff_time, diff_series)
plt.show()

# + id="B29KtYwRmBJU" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 388} outputId="755fa75d-4cbd-42b1-ebca-b40461368e2b"
# Use Moving Average function again after taking the difference in difference which is stationary.
diff_moving_avg = moving_average_forecast(diff_series, 50)[split_time - 365 - 50:]
plt.figure(figsize=(10, 6))
plot_series(time_valid, diff_series[split_time - 365:])
plot_series(time_valid, diff_moving_avg)
plt.show()

# + id="yNZIC2hymghc" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 388} outputId="b398d90d-b028-48ff-810f-f141fe02f375"
# Bring back trend and seasonality by adding the past values from 365 days ago.
diff_moving_avg_plus_past = series[split_time - 365:-365] + diff_moving_avg
plt.figure(figsize=(10, 6))
plot_series(time_valid, x_valid)
plot_series(time_valid, diff_moving_avg_plus_past)
plt.show()

# + id="RBGpyaO6mrla" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 52} outputId="126878c9-562c-437c-dd5f-396ebe9496e5"
# Compute MSE and MAE between (Moving Average plus past values) forecasts and predictions in the validation period. 
# This Moving Average is window size 30. 
print(keras.metrics.mean_squared_error(x_valid, diff_moving_avg_plus_past).numpy())
print(keras.metrics.mean_absolute_error(x_valid, diff_moving_avg_plus_past).numpy())

# + id="ofwTPbeBm_ey" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 388} outputId="09d6e177-eab2-4fdf-a7a6-47f16688c22e"
# Also use Moving Average of 10 windows on past values.
diff_moving_avg_plus_smooth_past =  moving_average_forecast(series[split_time - 370:-360], 10) + diff_moving_avg
plt.figure(figsize=(10, 6))
plot_series(time_valid, x_valid)
plot_series(time_valid, diff_moving_avg_plus_smooth_past)
plt.show()

# + id="X0tIaJAMm_KU" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 52} outputId="c3661fa6-0444-447e-f2b5-c26047520708"
# Compute MSE and MAE between (Moving Average on past values plus past values) forecasts and predictions in the validation period. 
# This Moving Average on past values is window size 10. 
print(keras.metrics.mean_squared_error(x_valid, diff_moving_avg_plus_smooth_past).numpy())
print(keras.metrics.mean_absolute_error(x_valid, diff_moving_avg_plus_smooth_past).numpy())


# + id="LJwMFE__YWUA" colab_type="code" colab={}
# 1st definition of autocorrelation.
def autocorrelation(time, amplitude, seed=None):
    rnd = np.random.RandomState(seed)
    φ1 = 0.5
    φ2 = -0.1
    ar = rnd.randn(len(time) + 50)
    ar[:50] = 100
    for step in range(50, len(time) + 50):
        ar[step] += φ1 * ar[step - 50]
        ar[step] += φ2 * ar[step - 33]
    return ar[50:] * amplitude


# + id="FiCUviRNYX9V" colab_type="code" colab={}
# 2nd definition of autocorrelation.
def autocorrelation(time, amplitude, seed=None):
    rnd = np.random.RandomState(seed)
    φ = 0.8
    ar = rnd.randn(len(time) + 1)
    for step in range(1, len(time) + 1):
        ar[step] += φ * ar[step - 1]
    return ar[1:] * amplitude


# + id="iBg-frAZZPJB" colab_type="code" colab={}
def impulses(time, num_impulses, amplitude=1, seed=None):
    rnd = np.random.RandomState(seed)
    impulse_indices = rnd.randint(len(time), size=10)
    series = np.zeros(len(time))
    for index in impulse_indices:
        series[index] += rnd.rand() * amplitude
    return series    


# + id="p7E8hcHkaH-4" colab_type="code" colab={}
def autocorrelation(source, φs):
    ar = source.copy()
    max_lag = len(φs)
    for step, value in enumerate(source):
        for lag, φ in φs.items():
            if step - lag > 0:
              ar[step] += φ * ar[step - lag]
    return ar


# + id="XtuVCLYwaKM0" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 283} outputId="c987f162-a3f9-4077-83e6-3986bac78887"
signal = impulses(time, 10, seed=42)
series = autocorrelation(signal, {1: 0.99})
plot_series(time, series)
plt.plot(time, signal, "k-")
plt.show()

# + id="nEPcqMgcaKtk" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 279} outputId="0e994cca-4a89-478f-b960-a5e8ade94afb"
signal = impulses(time, 10, seed=42)
series = autocorrelation(signal, {1: 0.70, 50: 0.2}) # The decay autocorrelates.
plot_series(time, series)
plt.plot(time, signal, "k-")
plt.show()

# + id="mKuMc3nabW5n" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 279} outputId="d9781c52-4a0b-450f-e2f4-7df1bc886550"
series_diff1 = series[1:] - series[:-1]
plot_series(time[1:], series_diff1)

# + id="BR0pmmgqbYxZ" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 301} outputId="93036f31-71a6-4c28-8adf-44f8a3e221ee"
from pandas.plotting import autocorrelation_plot

autocorrelation_plot(series)

# + id="VoGh6KYgbZbE" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 557} outputId="b9917b69-124d-45e6-d204-df5b289c8646"
from statsmodels.tsa.arima_model import ARIMA

model = ARIMA(series, order=(5, 1, 0))
model_fit = model.fit(disp=0)
print(model_fit.summary())

# + id="mZVby0dvcQVR" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 301} outputId="e67768be-d8d7-4e91-947e-bb23a71f7039"
series_diff = series
for lag in range(50):
  series_diff = series_diff[1:] - series_diff[:-1]

autocorrelation_plot(series_diff)

# + id="5b5EVPZpdBXz" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 265} outputId="2c9d7bfc-54bf-468f-f245-7d2044141add"
import pandas as pd
series_diff1 = pd.Series(series[1:] - series[:-1])
autocorrs = [series_diff1.autocorr(lag) for lag in range(1, 60)]
plt.plot(autocorrs)
plt.show()

# + id="4F1bQ0o_ekTZ" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 388} outputId="029c1001-f328-4f98-81d2-09ae6819bd16"
naive_forecast = series[split_time - 1:-1] # naive prediction is to assume value at t+1 is the same as t which is the current value

plt.figure(figsize=(10, 6))
plot_series(time_valid, x_valid)
plot_series(time_valid, naive_forecast)

# + id="-_Q1V7eGfAU5" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 388} outputId="33d0d174-9be4-4f1f-bf67-202ce4364985"
plt.figure(figsize=(10, 6))
plot_series(time_valid, x_valid, start=0, end=150)
plot_series(time_valid, naive_forecast, start=1, end=151)

StarcoderdataJupyterScriptsDedupFiltered