<jupyter_start><jupyter_text>Deep filtering<jupyter_code>import os
os.environ["CUDA_VISIBLE_DEVICES"] = "-1"
import collections

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

import keras
import tensorflow as tf
from tensorflow.python.keras import backend
# config = tf.compat.v1.ConfigProto(device_count={'GPU': 1, 'CPU': 6})
config = tf.compat.v1.ConfigProto(device_count={'CPU': 6})
sess = tf.compat.v1.Session(config=config)
backend.set_session(sess)<jupyter_output><empty_output><jupyter_text>We consider the basic econometric stochastic volatility model written in the form of a nonlinear state-space model: * <NAME>. (1990). *Forecasting, Structural Time Series Models and the Kalman Filter.* Cambridge University Press, New York. This model is no longer state-of-the-art, but we choose to start with something basic. First, we need some way to sample the parameters. We use the same priors as in * <NAME>., <NAME>. (2000). *BUGS for a Bayesian Analysis of Stochastic Volatility Models*.<jupyter_code>parameter_count = 3
def generate_parameters():
    precision = .1
    variance = 1. / precision
    std = np.sqrt(variance)
    mu = np.random.normal(loc=0., scale=std)
    phistar = np.random.beta(a=20., b=1.5)
    itau2 = np.random.gamma(shape=2.5, scale=.025)
    beta = np.exp(.5 * mu)
    phi = 2. * phistar - 1.
    tau = np.sqrt(1. / itau2)
    Parameters = collections.namedtuple('Parameters', ['mu', 'phi', 'tau'])
    return Parameters(mu=mu, phi=phi, tau=tau)<jupyter_output><empty_output><jupyter_text>Given the parameters, we can generate the latent states $x_0, x_1, \ldots, x_n$:<jupyter_code>def generate_xs(parameters, n):
    xs = []
    x0 = np.random.normal(loc=parameters.mu, scale=parameters.tau)
    xs.append(x0)
    for x in range(n):
        x_mean = parameters.mu + parameters.phi * (xs[-1] - parameters.mu)
        x = np.random.normal(loc=x_mean, scale=parameters.tau)
        xs.append(x)
    return np.array(xs)<jupyter_output><empty_output><jupyter_text>Given the latent states $x_0, x_1, \ldots, x_n$, we can generate the observations $y_1, y_2, \ldots, y_n$:<jupyter_code>def generate_ys(xs):
    return np.array([np.random.normal(loc=0., scale=np.exp(.5 * x)) for x in xs[1:]])<jupyter_output><empty_output><jupyter_text>Thus we have the machinery for generating observations using a simple Monte Carlo procedure. We fix<jupyter_code>n = 1000<jupyter_output><empty_output><jupyter_text>The neural networks that we will apply will operate on the basis of a sliding window of a given size, which we also fix<jupyter_code>window_size = 50
iter_count = 1000<jupyter_output><empty_output><jupyter_text>Before we start sampling random combinations of parameters, we fix them as in (Meyer, 2000):<jupyter_code>Parameters = collections.namedtuple('Parameters', ['mu', 'phi', 'tau'])
parameters = Parameters(mu=2. * np.log(.7204), phi=.9807, tau=.1489)
parameters<jupyter_output><empty_output><jupyter_text>Let's generate the latent states using these parameters,<jupyter_code>xs = generate_xs(parameters, n)
plt.plot(xs);<jupyter_output><empty_output><jupyter_text>Given the latent states, we can generate the observations:<jupyter_code>ys = generate_ys(xs)
plt.plot(ys);<jupyter_output><empty_output><jupyter_text>Our first question is whether, with the parameters fixed, a neural network can replace a particle filter. We first consider the architecture and procedure described in * <NAME>, <NAME>, <NAME> (2020). Deep Filtering. https://arxiv.org/abs/2008.03878<jupyter_code>model = keras.Sequential()
model.add(keras.layers.Dense(10, activation='sigmoid', input_dim=window_size))
model.add(keras.layers.Dense(10, activation='sigmoid'))
model.add(keras.layers.Dense(10, activation='sigmoid'))
model.add(keras.layers.Dense(10, activation='sigmoid'))
model.add(keras.layers.Dense(10, activation='sigmoid'))
model.add(keras.layers.Dense(1, activation='linear'))
keras.utils.plot_model(model, show_layer_names=False, show_shapes=True)<jupyter_output><empty_output><jupyter_text>However, we face a problem with this configuration. Let us first check whether this configuration can learn (effectively) the identity function in a trivial filtering problem$$x_{k+1} = x_k + 1,$$$$y_k = x_k.$$<jupyter_code>optimizer = keras.optimizers.SGD(lr=0.01)
model.compile(loss='mse', optimizer=optimizer, metrics=['mse'])
losses = []
xs = np.array(range(0, n+1))
yx = xs[1:]
xs = (xs - np.min(xs)) / (np.max(xs) - np.min(xs)) - .5
ys = (ys - np.min(ys)) / (np.max(ys) - np.min(ys)) - .5
for i in range(iter_count):
    inputs = []
    outputs = []
    for j in range(n - window_size):
        an_input = ys[j:j+window_size]
        inputs.append(an_input)
        an_output = xs[j+window_size]
        outputs.append(an_output)
    the_inputs = np.array(inputs)
    history = model.fit(the_inputs, np.array(outputs))
    losses.append(history.history['loss'][0])<jupyter_output>30/30 [==============================] - 1s 807us/step - loss: 0.8843 - mse: 0.8843
30/30 [==============================] - 0s 785us/step - loss: 0.0793 - mse: 0.0793
30/30 [==============================] - 0s 739us/step - loss: 0.0754 - mse: 0.0754
30/30 [==============================] - 0s 734us/step - loss: 0.0752 - mse: 0.0752
30/30 [==============================] - 0s 738us/step - loss: 0.0753 - mse: 0.0753
30/30 [==============================] - 0s 750us/step - loss: 0.0753 - mse: 0.0753
30/30 [==============================] - 0s 723us/step - loss: 0.0753 - mse: 0.0753
30/30 [==============================] - 0s 716us/step - loss: 0.0753 - mse: 0.0753
30/30 [==============================] - 0s 713us/step - loss: 0.0754 - mse: 0.0754
30/30 [==============================] - 0s 709us/step - loss: 0.0754 - mse: 0.0754
30/30 [==============================] - 0s 738us/step - loss: 0.0754 - mse: 0.0754
30/30 [==============================] - 0s 796us/step - loss: 0.0754 - mse:[...]<jupyter_text>Let us check the (in-sample) performance of the above learning procedure:<jupyter_code>plt.plot(losses)
plt.plot(pd.Series(losses).ewm(alpha=.1).mean());
pred_outputs = []
true_outputs = []
for j in range(n - window_size):
    an_input = np.array([ys[j:j+window_size]])
    an_output = xs[j+window_size]
    true_outputs.append(an_output)
    pred_output = model.predict(an_input)[0]
    pred_outputs.append(pred_output)
plt.plot(true_outputs, label='true_outputs')
plt.plot(pred_outputs, label='pred_outputs')
plt.legend();<jupyter_output><empty_output><jupyter_text>The reason for this problem is known. It is described and analysed in detail in * <NAME>, Hardt, Mozer, Singer (2020). Identity Crisis: Memorization and Generalization under Extreme Overparameterization. "CNNs of up to 10 layers are able to generalize from a single example, whereas FCNs cannot learn the identity function reliably from 60k examples." Let's, then, try a configuration based on CNNs.<jupyter_code>model = keras.Sequential()
model.add(keras.layers.Conv1D(filters=64, kernel_size=2, activation='relu', input_shape=(window_size, 1)))
model.add(keras.layers.Conv1D(filters=64, kernel_size=2, activation='relu'))
model.add(keras.layers.Conv1D(filters=64, kernel_size=2, activation='relu'))
model.add(keras.layers.Conv1D(filters=64, kernel_size=2, activation='relu'))
model.add(keras.layers.Flatten())
model.add(keras.layers.Dense(1, activation='linear'))
keras.utils.plot_model(model, show_layer_names=False, show_shapes=True)
optimizer = keras.optimizers.SGD(lr=0.01)
model.compile(loss='mse', optimizer=optimizer, metrics=['mse'])
model.save_weights('deep-filter-model.h5')
losses = []
for i in range(1000):
    xs = np.array(range(0, n+1))
    ys = xs[1:]
    xs = (xs - np.min(xs)) / (np.max(xs) - np.min(xs)) - .5
    ys = (ys - np.min(ys)) / (np.max(ys) - np.min(ys)) - .5
    inputs = []
    outputs = []
    for j in range(n - window_size):
        an_input = ys[j:j+window_size]
        inputs.append(an_input)
        an_output = xs[j+window_size]
        outputs.append(an_output)
    the_inputs = np.array(inputs)
    history = model.fit(the_inputs.reshape(the_inputs.shape[0], the_inputs.shape[1], 1), np.array(outputs))
    losses.append(history.history['loss'][0])
plt.plot(losses)
plt.plot(pd.Series(losses).ewm(alpha=.1).mean());
pred_outputs = []
true_outputs = []
for j in range(n - window_size):
    an_input = np.array([ys[j:j+window_size]])
    an_output = xs[j+window_size]
    true_outputs.append(an_output)
    pred_output = model.predict(np.reshape(an_input, (an_input.shape[0], an_input.shape[1], 1)))[0]
    pred_outputs.append(pred_output)
plt.plot(true_outputs, label='true_outputs')
plt.plot(pred_outputs, label='pred_outputs')
plt.legend();<jupyter_output><empty_output><jupyter_text>Much better! Now let's try this on the stochastic volatility model (with parameters fixed as above). We train the filter, as in (Zhang, 2020), using a MC simulation...<jupyter_code>model.load_weights('deep-filter-model.h5')
model.compile(loss='mse', optimizer=optimizer, metrics=['mse'])
losses = []
for i in range(1000):
    xs = generate_xs(parameters, n)
    ys = generate_ys(xs)
    xs = (xs - np.min(xs)) / (np.max(xs) - np.min(xs)) - .5
    ys = (ys - np.min(ys)) / (np.max(ys) - np.min(ys)) - .5
    inputs = []
    outputs = []
    for j in range(n - window_size):
        an_input = ys[j:j+window_size]
        inputs.append(an_input)
        an_output = xs[j+window_size]
        outputs.append(an_output)
    the_inputs = np.array(inputs)
    history = model.fit(the_inputs.reshape(the_inputs.shape[0], the_inputs.shape[1], 1), np.array(outputs))
    losses.append(history.history['loss'][0])
plt.plot(losses)
plt.plot(pd.Series(losses).ewm(alpha=.1).mean());<jupyter_output><empty_output><jupyter_text>How does the model perform out-of-sample?<jupyter_code>xs = generate_xs(parameters, n)
ys = generate_ys(xs)
xs = (xs - np.min(xs)) / (np.max(xs) - np.min(xs)) - .5
ys = (ys - np.min(ys)) / (np.max(ys) - np.min(ys)) - .5
pred_outputs = []
true_outputs = []
for j in range(n - window_size):
    an_input = np.array([ys[j:j+window_size]])
    an_output = xs[j+window_size]
    true_outputs.append(an_output)
    pred_output = model.predict(np.reshape(an_input, (an_input.shape[0], an_input.shape[1], 1)))[0]
    pred_outputs.append(pred_output)
plt.plot(true_outputs, label='true_outputs')
plt.plot(pred_outputs, label='pred_outputs')
plt.legend();<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1807806655797449, 'fraction_numerical': 0.03749880118922029, 'mean_word_length': 4.925, 'pattern_counts': {'":': 0, '<': 67, '<?xml version=': 0, '>': 67, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '317756', 'license': 'Apache-2.0', 'n_tokens_mistral': 4000, 'n_tokens_neox': 3673, 'n_words': 807, 'path': 'jupyter-notebooks/03-deep-filtering.ipynb', 'repo_name': 'sydx/deep-econometrics', 'size': '10448'}