In this module, we prepare the time series model for the uncertain inflow energy $X_t$.

Whole model description:  
https://www2.isye.gatech.edu/people/faculty/Alex_Shapiro/Rep-Dec2011.pdf


```python
import pandas
import numpy
# read historical data for the four reserviors
hist = [
    pandas.read_csv(
        "./data/hist_{}.csv".format(i), 
        sep=";",
        index_col=0,
    ) 
    for i in range(4)
]
# concatenate the data to fill missing value yearly
hist_all = pandas.concat(hist, axis=1, ignore_index=True)
hist_all.fillna(hist_all.mean(), inplace=True)
hist_all_log = numpy.log(hist_all)
# split the data back to individual reserviors
hist = [hist_all.iloc[:,12*i:12*(i+1)] for i in range(4)]
```

Let $Y_t = log(X_t)$. We hope to fit into the following model:  
\begin{equation}
    Y_t - \mu_t = \phi_t(Y_{t-1}-\mu_{t-1})+ \epsilon_t
\end{equation}
Where: 
1. $\mu_t=\mu_{t+12}$ be the monthly averages of $Y_t$  
2. $\epsilon_t$ independently follows multivariate normal distribution

The above model is not linear w.r.t $X_t$. Plug in $Y_t = log(X_t)$ and linearize the model (first order taylor expansion approximation), we get the following model  

\begin{equation}
    X_t = e^{\epsilon_t} [\tilde{\mu}_t + \gamma_t \frac{\tilde{\mu}_t}{\tilde{\mu}_{t-1}} (X_{t-1}-\tilde{\mu}_{t-1})]
\end{equation}
Where:
1. $\gamma_t = \gamma_{t+12}$  
2. $\tilde{\mu}_t = e^{\mu_t},\tilde{\mu}_{t+12}=\tilde{\mu}_t$    
2. $\epsilon_t$ independently follows multivariate normal distribution $N(0,\Sigma_t)$, $\Sigma_{t+12}=\Sigma_{t}$  

This model is called auto-regressive process with multiplicative error. We will fit into this model.  

Denote $R_t = \frac{X_t - \tilde{\mu}_t}{\tilde{\mu}_t}$ and set $\epsilon_t$ to 0, the above model can be written as  
\begin{equation}
    R_t = \gamma_t R_{t-1}
\end{equation}
$\Sigma_t$ is modelled as sample covariance matrix for  
\begin{equation}
    log(\frac{X_t}{(1-\gamma_t)\tilde{\mu}_t +\gamma_t \frac{\tilde{\mu}_t}{\tilde{\mu}_{t-1}} X_{t-1}})
\end{equation}

Compute $\tilde{\mu}_t$ 


```python
hist_log = [numpy.log(hist[i]) for i in range(4)]
mu = pandas.DataFrame([hist_log[i].mean().values for i in range(4)]).T
exp_mu = numpy.exp(mu)
exp_mu.head()
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>0</th>
      <th>1</th>
      <th>2</th>
      <th>3</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>54330.003557</td>
      <td>6255.342757</td>
      <td>13512.731607</td>
      <td>9902.926930</td>
    </tr>
    <tr>
      <th>1</th>
      <td>56588.831013</td>
      <td>7026.639444</td>
      <td>13582.911727</td>
      <td>13242.238584</td>
    </tr>
    <tr>
      <th>2</th>
      <td>53120.586175</td>
      <td>6086.382413</td>
      <td>13227.368691</td>
      <td>15756.100521</td>
    </tr>
    <tr>
      <th>3</th>
      <td>40532.704941</td>
      <td>5426.700993</td>
      <td>10897.449454</td>
      <td>15786.030067</td>
    </tr>
    <tr>
      <th>4</th>
      <td>29406.129246</td>
      <td>6187.208315</td>
      <td>6619.769923</td>
      <td>10521.234833</td>
    </tr>
  </tbody>
</table>
</div>



Perform least square fit to the $R_t$ sequence on a monthly basis. Beware to keep temporal order (December data should be lagged properly).


```python
gamma = {i:[None]*12 for i in range(4)}
for i in range(4):
    for t in range(12):
        if t == 0:
            X = hist[i].iloc[:-1,11]
            Y = hist[i].iloc[1:,0]
        else:
            X = hist[i].iloc[:,t-1]
            Y = hist[i].iloc[:,t]
        X = X/exp_mu.iloc[t-1,i] - 1
        Y = Y/exp_mu.iloc[t,i] - 1
        gamma[i][t] = numpy.dot(X,Y)/numpy.dot(X,X)
gamma = pandas.DataFrame.from_dict(gamma)
gamma.head()
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>0</th>
      <th>1</th>
      <th>2</th>
      <th>3</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>0.645308</td>
      <td>0.424965</td>
      <td>0.520670</td>
      <td>0.695171</td>
    </tr>
    <tr>
      <th>1</th>
      <td>0.589177</td>
      <td>0.659803</td>
      <td>0.844489</td>
      <td>0.560689</td>
    </tr>
    <tr>
      <th>2</th>
      <td>0.589719</td>
      <td>0.542662</td>
      <td>0.970522</td>
      <td>0.629378</td>
    </tr>
    <tr>
      <th>3</th>
      <td>0.712945</td>
      <td>0.848216</td>
      <td>0.588331</td>
      <td>0.670044</td>
    </tr>
    <tr>
      <th>4</th>
      <td>0.744759</td>
      <td>0.700799</td>
      <td>0.904786</td>
      <td>0.930524</td>
    </tr>
  </tbody>
</table>
</div>



Compute the sample cov


```python
sigma = [None] * 12
for t in range(12):
    if t == 0:
        X = pandas.DataFrame([hist[i].iloc[:-1,11] for i in range(4)]).T
        Y = pandas.DataFrame([hist[i].iloc[1:,0] for i in range(4)]).T
    else:
        X = pandas.DataFrame([hist[i].iloc[:,t-1] for i in range(4)]).T
        Y = pandas.DataFrame([hist[i].iloc[:,t] for i in range(4)]).T
    X.reset_index(drop=True,inplace=True)
    Y.reset_index(drop=True,inplace=True)
    X.columns = [0,1,2,3]
    Y.columns = [0,1,2,3]
    matrix = numpy.log(
        Y
        / (
            (1-gamma.iloc[t]) * exp_mu.iloc[t] 
            + gamma.iloc[t] * exp_mu.iloc[t]/exp_mu.iloc[t-1] * X
        )
    )
    sigma[t] = pandas.DataFrame(numpy.cov(matrix.T,ddof=1))
sigma[11]
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>0</th>
      <th>1</th>
      <th>2</th>
      <th>3</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>0.037904</td>
      <td>-0.002144</td>
      <td>0.020849</td>
      <td>0.014207</td>
    </tr>
    <tr>
      <th>1</th>
      <td>-0.002144</td>
      <td>0.222603</td>
      <td>-0.054403</td>
      <td>-0.060221</td>
    </tr>
    <tr>
      <th>2</th>
      <td>0.020849</td>
      <td>-0.054403</td>
      <td>0.088267</td>
      <td>0.056809</td>
    </tr>
    <tr>
      <th>3</th>
      <td>0.014207</td>
      <td>-0.060221</td>
      <td>0.056809</td>
      <td>0.073611</td>
    </tr>
  </tbody>
</table>
</div>



Now we can construct a sample path generator for the above time series model


```python
def generator(random_state):
    initial = [[41248.715300,7386.860854,10124.561460,6123.808537]]
    for t in range(12):
        noise = numpy.exp(random_state.multivariate_normal(mean=[0]*4, cov=sigma[t]))
        new = pandas.Series(noise) * (
            (1-gamma.iloc[t]) * exp_mu.iloc[t]
            + gamma.iloc[t] * exp_mu.iloc[t]/exp_mu.iloc[t-1] * initial[-1]
        )
        initial.append(new)
    del initial[0]    
    return initial
```

Use our sample path generator to simulate inflow data. 


```python
import matplotlib.pyplot as plt
from MSP.utils.plot import fan_plot
simulation = [[] for _ in range(4)]
for iter in range(100):
    s = generator(numpy.random.RandomState(iter))  
    for i in range(4):
        s0 = [item[i] for item in s]
        simulation[i].append(s0)
df_sim = [pandas.DataFrame(simulation[i]) for i in range(4)]
```

Use fan plots to compare simulated data and historical data. The first column shows the simulated data while the second column shows the historical data. It is clear that the shapes are similar, except that the simulated data for the second reservior have a bit more tails than the real data. This behavior is also recognized in the paper:  

http://www.optimization-online.org/DB_FILE/2017/09/6201.pdf  

Overall, the above time series model fits well to the data. 

For more discussion of validating the above model, have a look at the following paper:    
https://www2.isye.gatech.edu/people/faculty/Alex_Shapiro/Rep-Dec2011.pdf


```python
fig = plt.figure(figsize=(10,10))
ax = [None] * 4
ax[0] = plt.subplot(421)
ax[1] = plt.subplot(423)
ax[2] = plt.subplot(425)
ax[3] = plt.subplot(427)
for i in range(4):
    fan_plot(df_sim[i], ax[i])
bx = [None] * 4
bx[0] = plt.subplot(422, sharey=ax[0])
bx[1] = plt.subplot(424, sharey=ax[1])
bx[2] = plt.subplot(426, sharey=ax[2])
bx[3] = plt.subplot(428, sharey=ax[3])
for i in range(4):
    fan_plot(hist[i], bx[i])
```

Output coefficients of the generator for later use.


```python
gamma.to_csv("./data/gamma.csv")
for t in range(12):
    sigma[t].to_csv("./data/sigma_{}.csv".format(t))
exp_mu.to_csv("./data/exp_mu.csv")
```
