text:
```python
# Perform PCA, classify w/ log regression+ovr, figure decision boundaries
# Jan Nagler (adapted)
# applied to wine database
# always works
import pandas as pd
from sklearn.datasets import make_blobs
from numpy import mean
from numpy import std
from numpy.random import randn
from numpy.random import seed
from numpy import cov
from scipy.stats import pearsonr
from matplotlib import pyplot
from matplotlib.colors import ListedColormap
import sympy as sympy
import numpy as np
from scipy import *
import matplotlib.pyplot as plt
from sklearn.cluster import KMeans
import seaborn as sns
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from IPython.display import display, Math, Latex
from sklearn.linear_model import LogisticRegression
from sklearn.decomposition import PCA
```

```python
# Load a database for wines, labled (3 classes), 13 attributes, N=178
# See www page: https://archive.ics.uci.edu/ml/datasets/wine for description of data
#1) Alcohol
#2) Malic acid
#3) Ash
#4) Alcalinity of ash
#5) Magnesium
#6) Total phenols
#7) Flavanoids
#8) Nonflavanoid phenols
#9) Proanthocyanins
#10)Color intensity
#11)Hue
#12)OD280/OD315 of diluted wines
#13)Proline

# Create data frame (data + labels) from www source
df_wine = pd.read_csv('https://archive.ics.uci.edu/ml/'
                      'machine-learning-databases/wine/wine.data',
                      header=None)

# Display dimensions of dataset
print( df_wine.shape )

# Display first rows (head) of dataset
df_wine.head(5)
```

    (178, 14)

      0
      1
      2
      3
      4
      5
      6
      7
      8
      9
      10
      11
      12
      13

      0
      1
      14.23
      1.71
      2.43
      15.6
      127
      2.80
      3.06
      0.28
      2.29
      5.64
      1.04
      3.92
      1065

      1
      1
      13.20
      1.78
      2.14
      11.2
      100
      2.65
      2.76
      0.26
      1.28
      4.38
      1.05
      3.40
      1050

      2
      1
      13.16
      2.36
      2.67
      18.6
      101
      2.80
      3.24
      0.30
      2.81
      5.68
      1.03
      3.17
      1185

      3
      1
      14.37
      1.95
      2.50
      16.8
      113
      3.85
      3.49
      0.24
      2.18
      7.80
      0.86
      3.45
      1480

      4
      1
      13.24
      2.59
      2.87
      21.0
      118
      2.80
      2.69
      0.39
      1.82
      4.32
      1.04
      2.93
      735

```python
# Read out samples (X) and labels (y)
X, y = df_wine.iloc[:, 1:].values, df_wine.iloc[:, 0].values

# Split into training and testing sets
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.3, stratify=y, random_state=0)

# Standardize the features (zero mean, unit variance)
sc = StandardScaler()

# Fit results must be used later (mu and sigma)
X_train_std = sc.fit_transform(X_train)

# Normalize test data set with mu/sigma of training data
X_test_std = sc.transform(X_test)
```

```python
# Compute covariance matrix and eigenvalues (EVal) / eigenvectors (EVec)
print (X_train.shape)
cov_mat = np.cov(X_train_std.T) #cov matrix from data
EVal, EVec = np.linalg.eig(cov_mat)
```

    (124, 13)

```python
# Calculate (sum of) explained variances
sum_EVal = sum(EVal)
var_exp = [(i / sum_EVal) for i in sorted(EVal, reverse=True)]

# Plot explained variances
plt.bar(range(1,14), var_exp, alpha=0.3,
        align='center', label='individual explained variance')
plt.ylabel('Explained variance ratio')
plt.xlabel('Principal component index')
plt.legend(loc='best')
plt.show()
```

```python
# Make a list of (eigenvalue, eigenvector) tuples
eigen_pairs = [(np.abs(EVal[i]), EVec[:, i]) for i in range(len(EVal))]

# Sort the (EVal, EVec) tuples from high to low (reverse), key=EV
eigen_pairs.sort(key=lambda k: k[0], reverse=True)
```

```python
# Construct matrix w
# stacking, reshaping, adding dim of arrays
w = np.hstack((eigen_pairs[0][1][:, np.newaxis], eigen_pairs[1][1][:, np.newaxis]))
print(w)
#the above code is equivalent to that one (using reshape)
#w = np.hstack((eigen_pairs[0][1].reshape(13,1), eigen_pairs[1][1].reshape(13,1)))
#print(w)
```

    [[-0.13724218  0.50303478]
     [ 0.24724326  0.16487119]
     [-0.02545159  0.24456476]
     [ 0.20694508 -0.11352904]
     [-0.15436582  0.28974518]
     [-0.39376952  0.05080104]
     [-0.41735106 -0.02287338]
     [ 0.30572896  0.09048885]
     [-0.30668347  0.00835233]
     [ 0.07554066  0.54977581]
     [-0.32613263 -0.20716433]
     [-0.36861022 -0.24902536]
     [-0.29669651  0.38022942]]

```python
# Transform the entire (training) data to the PCA space, X'=X.W
X_train_pca = X_train_std.dot(w)
print (X_train_std.shape, w.shape)
```

    (124, 13) (13, 2)

```python
# Display LaTeX if you want (but nobody wants that)
display(Math(r'X^*=XW'))
```

$$X^*=XW$$

```python
# Visualize the data
colors = ['r', 'b', 'g']
markers = ['s', 'v', 'o']
#
for l, c, m in zip(np.unique(y_train), colors, markers):
    plt.scatter(X_train_pca[y_train==l, 0], # Pick x component for label
                X_train_pca[y_train==l, 1],
                color=c, label=l, marker=m, alpha=0.4)
plt.xlabel('PC 1')
plt.ylabel('PC 2')
plt.legend(loc='lower left')
plt.show()
```

```python
# Set up PCA and logistic regression model
pca = PCA(n_components=2)
lr = LogisticRegression(multi_class='ovr', solver='liblinear')

# Fit and transform training data, given on PCA reduction to k(=2) principle components
X_train_pca = pca.fit_transform(X_train_std)
X_test_pca = pca.transform(X_test_std)
# solves task, given 3 classes (as from y_train)
lr.fit(X_train_pca, y_train)
```

    LogisticRegression(C=1.0, class_weight=None, dual=False, fit_intercept=True,
                       intercept_scaling=1, l1_ratio=None, max_iter=100,
                       multi_class='ovr', n_jobs=None, penalty='l2',
                       random_state=None, solver='liblinear', tol=0.0001, verbose=0,
                       warm_start=False)

```python
def plot_decision_regions(X, y, classifier, resolution=0.01):
    # Does not work....always....Setup marker generator and color map: FIXED!
    markers = ('s', 'v', 'o')
    colors = ('r', 'b', 'g')
    cmap = ListedColormap(colors[:len(np.unique(y))])

    # plot the decision surface
    x1_min, x1_max = X[:, 0].min() - 1, X[:, 0].max() + 1
    x2_min, x2_max = X[:, 1].min() - 1, X[:, 1].max() + 1
    xx1, xx2 = np.meshgrid(np.arange(x1_min, x1_max, resolution),
                           np.arange(x2_min, x2_max, resolution))

    # Z is the prediction of the class, given point in plane
    Z = classifier.predict(np.array([xx1.ravel(), xx2.ravel()]).T)
    Z = Z.reshape(xx1.shape)

    # Z=f(xx1,yy1), plot classes in plane using color map but opaque
    plt.contourf(xx1, xx2, Z, alpha=0.2, cmap=cmap)
    plt.xlim(xx1.min(), xx1.max())
    plt.ylim(xx2.min(), xx2.max())

    # Plot data points, given labels
    for idx, cl in enumerate(np.unique(y)):
        plt.scatter(x=X[y == cl, 0],
                    y=X[y == cl, 1],
                    alpha=0.6,
                    c=[cmap(idx)],
                    edgecolor='black',
                    marker=markers[idx],
                    label=cl)
```

```python
plot_decision_regions(X_train_pca, y_train, classifier=lr)
plt.xlabel('PC 1')
plt.ylabel('PC 2')
plt.legend(loc='lower left')
plt.show()
```

```python
# plot decision regions for test set
plot_decision_regions(X_test_pca, y_test, classifier=lr)
plt.xlabel('PC1')
plt.ylabel('PC2')
plt.legend(loc='lower left')
plt.show()
```

```python
print("Hooray")
```

    Hooray

```python

```


text:
```python
from sympy.physics.units import *
from sympy import *

# Rounding:
import decimal
from decimal import Decimal as DX
def iso_round(obj, pv, rounding=decimal.ROUND_HALF_EVEN):
    import sympy
    """
    Rounding acc. to DIN EN ISO 80000-1:2013-08
    place value = Rundestellenwert
    """
    assert pv in set([
        # place value   #  round to:
        100,            #  3rd last digit before decimal
        10,             #  2nd last
        1,              #  last
        0.1,            #  1st digit after decimal
        0.01,           #  2nd
        0.001,          #  3rd
        0.0001,         #  4th
        0.00001,        #  5th
        0.000001,       #  6th
        0.0000001,      #  7th
        0.00000001,     #  8th
        0.000000001,    #  9th
        0.0000000001,   # 10th
        ])
    try:
        tmp = DX(str(float(obj)))
        obj = tmp.quantize(DX(str(pv)), rounding=rounding)
    except:
        for i in range(len(obj)):
            tmp = DX(str(float(obj[i])))
            obj[i] = tmp.quantize(DX(str(pv)), rounding=rounding)
    return obj

# LateX:
kwargs = {}
kwargs["mat_str"] = "bmatrix"
kwargs["mat_delim"] = ""
# kwargs["symbol_names"] = {FB: "F^{\mathsf B}", }

# Units:
(k, M, G ) = ( 10**3, 10**6, 10**9 )
(mm, cm, deg) = ( m/1000, m/100, pi/180)
Newton = kg*m/s**2
Pa     = Newton/m**2
MPa    = M*Pa
GPa    = G*Pa
kN     = k*Newton

# ---

c1, c2, m1, m2 = var("c1, c2, m1, m2", positive=True)

c, mass = var("c, mass", positive=True)

K = Matrix([
[c1 + c2,  -c2],
[-c2,  c2],
])

M = Matrix([
[m1, 0],
[0, m2],
])

sub_list = [
    (m1,mass),
    (m2,2*mass),
    (c1,c),
    (c2,c),
]

# ξ = λ²
xi = var("xi")
w = var("omega")

pprint("\nK + λ² M:")
A = K + xi*M
pprint(A)

pprint("\nK + λ²M for (m1=m, m2=2m, c1=c, c2=c):")
A = A.subs(sub_list)
pprint(A)
# exit()
pprint("\nCharacteristic equation:")
eq = Eq(det(A))
pprint(eq)

sol_xi = solve(eq,xi)

# Init:
eigenfreqs = [0,0]
eigenvects = [Matrix([1,0]), Matrix([1,0])]

x1, x2 = var("x1, x2")
w = Matrix([x1, x2])

pprint("\n")
for i in range(len(sol_xi)):
    pprint("\n--- (Value, Vector)-Pair Number "+str(i+1)+" ---")
    pprint("\nValue λ² / (c/m):")
    xii = sol_xi[i]
    eigenfreqs[i] = sqrt(-xii)
    tmp = xii / (c/mass)
    pprint(tmp)
    Ai = A.subs(xi,xii)
    pprint("\nVector:")
    eq = Eq(Ai*w,Matrix([0,0]))
    sol = solve(eq, x2)
    pprint(sol)
    tmp = sol[x2]
    tmp = tmp.subs(x1, 1)
    eigenvects[i][1] = tmp

pprint("\nEigenfrequencies / sqrt(c/m) and resp. Eigenvector for x1=1:")
for i in range(2):
    pprint("\n--- (Value, Vector)-Pair Number "+str(i+1)+" ---")
    tmp = eigenfreqs[i] / sqrt(c/mass)
    pprint(tmp)
    tmp = iso_round(tmp,0.01)
    print("≈ ", end="")
    pprint(tmp)
    pprint("\n")
    pprint(eigenvects[i])

sub_list_ex = [
    (c,1*Newton/m),
    (mass, 1*kg),
    ]

# Set Me!
first = True  # first eigenmode
first = False # second eigenmode

if first:
    pprint("\nFirst Eigenmode:")
    tmp = eigenfreqs[0]
    fac = eigenvects[0][1]
else:
    pprint("\nSecond Eigenmode:")
    tmp = eigenfreqs[1]
    fac = eigenvects[1][1]

pprint("\nAngular Frequency / (1/s):")
tmp = tmp.subs(sub_list)
tmp = tmp.subs(sub_list_ex)
tmp /= (1/s)
w_ex = N(tmp,15)
tmp = iso_round(tmp,0.01)
pprint(tmp)
pprint("\nTime Period / s:")
tmp = 2*pi/w_ex
tmp = iso_round(tmp,0.01)
pprint(tmp)

from pylab import *
t_in_s = linspace(0,15,100)
cos_arg = w_ex  * t_in_s
x1_in_m = array([cos(float(x)) for x in cos_arg])
x2_in_m = fac * x1_in_m

tmp = plt.axis()
tmp = plt.grid()
tmp = plt.plot(t_in_s, x1_in_m, label=r"$x_1\,\, /  \,\, \mathrm{m}$")
tmp = plt.plot(t_in_s, x2_in_m, label=r"$x_2\,\, /  \,\, \mathrm{m}$")
tmp = plt.xlabel(r"$t  \,\, /  \,\, \mathrm{s}$")
tmp = plt.legend()
plt.show()

#
# K + λ² M:
# ⎡c₁ + c₂ + m₁⋅ξ     -c₂   ⎤
# ⎢                         ⎥
# ⎣     -c₂        c₂ + m₂⋅ξ⎦
#
# K + λ²M for (m1=m, m2=2m, c1=c, c2=c):
# ⎡2⋅c + mass⋅ξ       -c     ⎤
# ⎢                          ⎥
# ⎣     -c       c + 2⋅mass⋅ξ⎦
#
# Characteristic equation:
#    2
# - c  + (c + 2⋅mass⋅ξ)⋅(2⋅c + mass⋅ξ) = 0
#
#
# --- (Value, Vector)-Pair Number 1 ---
#
# Value λ² / (c/m):
#   5   √17
# - ─ + ───
#   4    4
#
# Vector:
# ⎧    x₁⋅(3 + √17)⎫
# ⎨x₂: ────────────⎬
# ⎩         4      ⎭
#
# --- (Value, Vector)-Pair Number 2 ---
#
# Value λ² / (c/m):
#   5   √17
# - ─ - ───
#   4    4
#
# Vector:
# ⎧    x₁⋅(-√17 + 3)⎫
# ⎨x₂: ─────────────⎬
# ⎩          4      ⎭
#
# Eigenfrequencies / sqrt(c/m) and resp. Eigenvector for x1=1:
#
# --- (Value, Vector)-Pair Number 1 ---
#   __________
# ╲╱ -√17 + 5
# ────────────
#      2
# ≈ 0.47
#
# ⎡   1   ⎤
# ⎢       ⎥
# ⎢3   √17⎥
# ⎢─ + ───⎥
# ⎣4    4 ⎦
# ≈
# ⎡1.0 ⎤
# ⎢    ⎥
# ⎣1.78⎦
#
# --- (Value, Vector)-Pair Number 2 ---
#   _________
# ╲╱ √17 + 5
# ───────────
#      2
# ≈ 1.51
#
# ⎡    1    ⎤
# ⎢         ⎥
# ⎢  √17   3⎥
# ⎢- ─── + ─⎥
# ⎣   4    4⎦
# ≈
# ⎡ 1.0 ⎤
# ⎢     ⎥
# ⎣-0.28⎦
#
# Example:
#
# Angular Frequency / (1/s):
# 0.47
#
# Time Period / s:
# 13.42

```

ext: ipynb

max_stars_repo_name: kassbohm/tm-snippets

max_issues_repo_name: kassbohm/tm-snippets

max_forks_repo_name: kassbohm/tm-snippets


text:
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()
```

      0
      1
      2
      3

      0
      54330.003557
      6255.342757
      13512.731607
      9902.926930

      1
      56588.831013
      7026.639444
      13582.911727
      13242.238584

      2
      53120.586175
      6086.382413
      13227.368691
      15756.100521

      3
      40532.704941
      5426.700993
      10897.449454
      15786.030067

      4
      29406.129246
      6187.208315
      6619.769923
      10521.234833

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()
```

      0
      1
      2
      3

      0
      0.645308
      0.424965
      0.520670
      0.695171

      1
      0.589177
      0.659803
      0.844489
      0.560689

      2
      0.589719
      0.542662
      0.970522
      0.629378

      3
      0.712945
      0.848216
      0.588331
      0.670044

      4
      0.744759
      0.700799
      0.904786
      0.930524

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]
```

      0
      1
      2
      3

      0
      0.037904
      -0.002144
      0.020849
      0.014207

      1
      -0.002144
      0.222603
      -0.054403
      -0.060221

      2
      0.020849
      -0.054403
      0.088267
      0.056809

      3
      0.014207
      -0.060221
      0.056809
      0.073611

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")
```

max_stars_repo_name: odow/msppy

max_issues_repo_name: odow/msppy

max_forks_repo_name: odow/msppy


text:
# Hopfield网络

Hopfield不同于BP网络，Hopfield网络具有反馈连接，即Hopfield应当属于递归神经网络。

Hopfield网络可由硬件实现。其可分为连续型和离散型两类，分别命名为CHNN(Continues Hopfield Neural Network)和DHNN(Discrete Hopfield Neural Network)。这里主要讨论的是离散型

## DHNN
离散型Hopfield网络是单层反馈网络，每一个节点的输出均会反馈至输入。

### 定义符号：

|符号|含义|
|:-:|:-:|
|$\bm{x(t)}$|t时刻的状态向量|
|$x_i$|第i个位置的状态值|
|$x_i^j$|第j个状态向量的i个位置的状态值|
|$\bm{W}$|连接权矩阵|
|$\bm{b}$|偏置|
|$f$|激活函数|
|$n$|状态维度|
|$k$|神经元数量|

考虑双极型离散数据，即状态只能是-1或1

### DHNN的基本概念
$$
\begin{equation}
            \bm{x(t+1)} = sign(\bm{W}x(t) + \bm{b})
\end{equation}
$$

对于Hopfield网络，首先需要指明初始状态$\bm{x(0)}$。然后按照上式进行反馈运行，网络最终达到某一种状态。对于状态有如下概念

1. **网络的稳定性**
若网络从初始状态开始，经过有限时间t后，网络的状态不再发生变化则称网络是稳定的
2. **网络的吸引子**
若网络达到稳定，则称稳定时的状态为网络的稳定吸引子；若网络状态有规律的在有限个状态之间振荡，则称网络处于有限环（limited circle）状态；若网络无规律的在某些状态之间振荡，则称网络处于混沌（chaos）状态
3. **吸引域**
对于某些特定的初始状态，网络按照一定的运行规则最后会稳定在同一个吸引子上。称能稳定在状态$x(t)$的所有初始状态集合为状态$x(t)$的吸引域

### 离散Hopfield网络的异步更新与全并行

* **异步更新**
对于异步更新方式，每一个时间步仅有某一个位置的状态更新，其余位置的状态不发生变化。对于这种更新方式有
$$
\begin{equation}
    x_i(t+1) = sign(\sum_{j=1}^n W_{ij} x_j(t) + b_i), i \in{1, 2, \cdots, n}
\end{equation}
$$

* **全并行工作方式**
全并行工作方式即式-1所示结果

### Hopfield网络的稳定性

* 能量函数
$$
\begin{equation}
    \begin{split}
        E
        &= -\frac{1}{2}\sum_i\sum_j W_{ij}x_ix_j + \sum_i b_i x_i \\
        &= -\frac{1}{2}\bm{x}^T \bm{W} \bm{x} + \bm{b}^T \bm{x}
    \end{split}
\end{equation}
$$

* **定理1**
若Hopfield网络为离散形式，且其状态按照异步的方式进行更新。假设连接权矩阵对称且对角线元素全为0（无自反馈），则网络状态会在有限步内收敛到稳定点

    证明略（以后再补吧）

* **定理2**
若Hopfield网络为连续形式，连接权矩阵对称且负定，则对于状态的每次同步变化，能量函数单调下降

    证明略（以后再补吧）

### Hopfield网络解决的问题

#### 联想记忆
若将稳态视为一个记忆，则由初态向稳态收敛的过程就是寻找记忆的过程。在这种观点下，初态是给定的部分信息，收敛过程即为从部分信息找到全部信息的过程（恢复过程），例如从噪声中恢复原有信息。

#### 优化计算
若将稳态视为某一优化问题目标函数的极小点，则由初态向稳态收敛的过程就是优化计算过程。网络的稳定状态即为对应问题的解

**对比：**
|视角|连接权矩阵$\bm{W}$|稳定状态|
|--|--|--|--|
|联想记忆|待学习|给定|
|优化计算|给定|待求解|

### Hopfield网络的参数学习

#### **联想记忆的设计方法**

对于DHNN，给定m个模式样本，且样本的维度为p：

$$
\begin{equation}
    \begin{array}{cc}
        S = [S^1, S^2, \cdots, S^m] \\
        \bm{\bm{S^k}} = (\bm{S^k}_1, \bm{S^k}_2, \cdots, \bm{S^k}_p)^T, k\in\{1, 2, \cdots, m\}
    \end{array}
\end{equation}
$$

若满足$S_i \in \{ -1, 1\}$并且模式样本两两正交可得连接权矩阵为
$$
\begin{equation}
    \bm{W} = \sum_{k=1}^m \bm{S^k} (\bm{S^k})^T - diag(\bm{S^k})^2
\end{equation}
$$

由于状态的取值只可能是-1, 1，上式可化简为
$$
\begin{equation}
    \bm{W} = \sum_{k=1}^m \bm{S^k} (\bm{S^k})^T - mI
\end{equation}
$$
其中$I$为单位矩阵

**证明**
对于上述确定的连接权矩阵，模式样本正好是网络的吸引子

由于模式样本两两正交，并且状态的取值只可能是-1, 1，因此有：
$$
\begin{equation}
    (\bm{S^k})^T\bm{S^l} =
    \left\{
        \begin{array}{cc}
        0&k\neq l \\
        ||\bm{S^k}||_2^2 = p &k=l
        \end{array}
    \right.
\end{equation}
$$

若初始状态为模式样本，则有:

$$
\begin{equation}
    x(0) = \bm{S^i}, i = 1, 2, \cdots, m
\end{equation}
$$

$$
\begin{equation}
    \begin{split}
        x(1)
        &= sign(\bm{W}x(0)) \\
        &= sign(\sum_{k=1}^m (\bm{S^k} (\bm{S^k})^T - mI) \bm{S^i}), i = 1, 2, \cdots, m \\
        &= sign(\sum_{k=1}^m(\bm{S^k} (\bm{S^k})^T \bm{S^i} - mI \bm{S^i})), i = 1, 2, \cdots, m \\
        &= sign(\bm{S^i} (\bm{S^i})^T \bm{S^i} - \sum_{k=1}^m(\bm{S^i})), i = 1, 2, \cdots, m \\
        &= sign((p-m)\bm{S^i}), i = 1, 2, \cdots, m
    \end{split}
\end{equation}
$$
由于模式样本两两正交，而$n$维空间最多有$n$个线性无关向量，因此样本的维度$p$必大于等于模式样本数量$m$，否则模式样本两两正交不成立，由此可得
$$
\begin{equation}
    x(1) = sign((p-m)\bm{S^i}) = sign(\bm{S^i}) = \bm{S^i} = x(0), i = 1, 2, \cdots, m
\end{equation}
$$

模式样本为吸引子

证毕

**注意！！！**
**上述证明仅证明了使用这种方式得到的权重可使得模式样本为吸引子，但是并没有保证吸引子仅有$m$个模式样本，即可能出现非模式样本的吸引子。非模式样本的吸引子也被称为“伪吸引子”**

```python
%matplotlib inline
import io
import os
import imageio
import cv2
from PIL import Image
import matplotlib.pyplot as plt
import matplotlib as mpl
import numpy as np
import random
```

```python
# hopfield
class Hopfield(object):

    def __init__(self, state_dims):

        # 维度
        self.dims = state_dims

        # 定义权重
        self._weight = np.zeros((state_dims, state_dims))

        # 记忆的模式
        self.restore_patterns = list()

        # 最大迭代次数
        self.max_steps = self.dims

    def store_pattern(self, pattern):
        self.restore_patterns.append(np.array(pattern))

    def weight_update(self):
        assert len(self.restore_patterns) > 0
        self._weight += np.multiply(self.restore_patterns[-1], self.restore_patterns[-1].T) - np.eye(self.dims)

    def _update_once(self, state_tensor, update_index, asynchronous=False):
        if asynchronous:
            temp_result = np.sum(self._weight[update_index, :] * state_tensor)
            if temp_result >= 0:
                new_value = 1
            else:
                new_value = -1
            if new_value != state_tensor[update_index]:
                change_value = True
            else:
                change_value = False
            state_tensor[update_index] = new_value
        else:
            temp_result = np.matmul(self._weight, state_tensor) >= 0
            temp_result = temp_result*2 - 1

            if np.sum(temp_result-state_tensor) != 0:
                change_value = True
            else:
                change_value = False

            state_tensor = temp_result
        return state_tensor, change_value

    def pred(self, state_tensor, asynchronous=False):
        temp_state_tensor = state_tensor.copy()
        steps = 0
        change_count = 0
        update_trace = list()
        update_trace.append(temp_state_tensor)

        while True:
            state_tensor, whether_change = self._update_once(state_tensor, steps%self.dims, asynchronous)

            if whether_change:
                change_count = 0
            else:
                change_count += 1

            if change_count == self.dims-200 or steps == self.max_steps:
                break

            temp_state_tensor = state_tensor
            update_trace.append(temp_state_tensor.copy())

            steps += 1

        return state_tensor, update_trace

def create_gif(output_name, image_num):
        """
        生成gif
        """
        with imageio.get_writer("./images/{}.gif".format(output_name), mode="I", fps=20) as Writer:
            for ind in range(image_num):
                image = imageio.imread("./images/temp/{}.png".format(ind))
                os.remove("./images/temp/{}.png".format(ind))
                Writer.append_data(image)
```

```python
# 测试
# 假设为状态只有两维
# 两个待记忆状态(-1, -1)和(1, 1)
# 测试输入只有四种情况(-1, -1), (-1, 1), (1, -1), (1, 1)

np.random.seed(1024)
random.seed(1024)

# --------------------------- 数据和模型 ----------------------------
# 待记忆模式
pattern_1 = np.array([-1, -1])
pattern_2 = np.array([1, 1])
# 测试输入
test_num = 4
test_data = np.array([[-1, -1], [-1, 1], [1, -1], [1, 1]])
# 定义模型
MyHopfield = Hopfield(2)

# --------------------------- 绘图相关 ------------------------------
# 迭代轨迹
update_traces = list()
# 最长周期
max_update_steps = 0
colors_list = list(mpl.colors.cnames.keys())
color_index_list = np.random.randint(0, len(colors_list), test_num)

# --------------------------- 正式流程 ------------------------------
# 记忆模式
MyHopfield.store_pattern(pattern_1)
MyHopfield.weight_update()
MyHopfield.store_pattern(pattern_2)
MyHopfield.weight_update()

# 联想记忆
for i in range(test_num):
    result_state, temp_trace = MyHopfield.pred(test_data[i].copy(), asynchronous=True)
    update_traces.append(temp_trace)

    print("状态{}对应的吸引子为{}".format(test_data[i], result_state))

    if len(update_traces[-1]) > max_update_steps:
        max_update_steps = len(update_traces[-1])

# 绘制
plt.figure(figsize=(8, 8), dpi=100)
for i in range(test_num):
    x_list = np.array(update_traces[i])[:, 0]
    y_list = np.array(update_traces[i])[:, 1]

    plt.scatter([x_list[0]], [y_list[0]], label="{}".format(test_data[i]))
    plt.plot(x_list, y_list)
    plt.legend()

plt.show()
```

```python
# 定义30*30的矩阵
# 记忆模式T、T_2
np.random.seed(1024)
random.seed(1024)

# --------------------------- 数据和模型 ----------------------------
# 待记忆模式
pattern_dims = 14

# 字母a
pattern_a = np.zeros((pattern_dims, pattern_dims, 1), dtype=np.uint8)
pattern_a = cv2.putText(pattern_a, "A", (0, 15), cv2.FONT_HERSHEY_SIMPLEX, 0.8, (255), 1)
pattern_a = np.array(pattern_a/127.5 - 1, dtype=np.int32).reshape((pattern_dims, pattern_dims))

# 字母b
pattern_b = np.zeros((pattern_dims, pattern_dims, 1), dtype=np.uint8)
pattern_b = cv2.putText(pattern_b, "B", (0, 13), cv2.FONT_HERSHEY_SIMPLEX, 0.6, (255), 1)
pattern_b = np.array(pattern_b/127.5 - 1, dtype=np.int32).reshape((pattern_dims, pattern_dims))

# 字母c
pattern_c = np.zeros((pattern_dims, pattern_dims, 1), dtype=np.uint8)
pattern_c = cv2.putText(pattern_c, "c", (0, 13), cv2.FONT_HERSHEY_SIMPLEX, 0.8, (255), 2)
pattern_c = np.array(pattern_c/127.5 - 1, dtype=np.int32).reshape((pattern_dims, pattern_dims))

print("模式1")
plt.imshow(pattern_a, cmap="gray")
plt.show()
print("模式2")
plt.imshow(pattern_b, cmap="gray")
plt.show()
print("模式3")
plt.imshow(pattern_c, cmap="gray")
plt.show()

# 测试输入, 前三个由字母T随机生成，后三个由字母T_2随机生成
test_num = 9
# 随机取值的像素数量
noise_num = 20
# 测试数据
test_data = list()
for i in range(test_num):
    random_position_x = np.random.randint(0, pattern_dims, noise_num)
    random_position_y = np.random.randint(0, pattern_dims, noise_num)

    if i <= 2:
        temp_test_matrix = pattern_a.copy()
    elif i <= 5:
        temp_test_matrix = pattern_b.copy()
    else:
        temp_test_matrix = pattern_c.copy()

    temp_test_matrix[random_position_x, random_position_y] =  - temp_test_matrix[random_position_x, random_position_y]
    test_data.append(temp_test_matrix)

# 定义模型
MyHopfield = Hopfield(pattern_dims**2)

# --------------------------- 绘图相关 ------------------------------
# 迭代轨迹
update_traces = list()

# --------------------------- 正式流程 ------------------------------
# 记忆模式
MyHopfield.store_pattern(pattern_a.flatten().reshape(-1, 1))
MyHopfield.weight_update()
MyHopfield.store_pattern(pattern_b.flatten().reshape(-1, 1))
MyHopfield.weight_update()
MyHopfield.store_pattern(pattern_c.flatten().reshape(-1, 1))
MyHopfield.weight_update()

print("权重矩阵")
plt.imshow(MyHopfield._weight, cmap="gray")
plt.show()

# 联想记忆
plt.figure(figsize=(10, 30), dpi=100)
for i in range(test_num):
    result_state, temp_trace = MyHopfield.pred(test_data[i].copy().flatten().reshape(-1, 1))
    update_traces.append(temp_trace)

    # 绘图
    plt.subplot(test_num, 2, 1+i*2)
    plt.title("input")
    plt.imshow(test_data[i].reshape(pattern_dims, pattern_dims), cmap="gray")
    plt.subplot(test_num, 2, 2+i*2)
    plt.title("output")
    plt.imshow(result_state.reshape(pattern_dims, pattern_dims), cmap="gray")

plt.show()

# # ----------------------------- 生成动图 --------------------------
# for i in range(test_num):
#     trace_ = update_traces[i]

#     print("create gif of {}-th input".format(i))
#     with imageio.get_writer("./images/hopfield_{}.gif".format(i), mode="I", fps=20) as Writer:

#         for count, step_img in enumerate(trace_):
#             plt.cla()
#             plt.clf()
#             fig = plt.figure("", frameon=False)
#             plt.title("step: {}".format(count))
#             plt.imshow(step_img.reshape(pattern_dims, pattern_dims), cmap="gray")
#             canvas = fig.canvas

#             # 创建输出流对象
#             buffer = io.BytesIO()
#             # 打印到输出流
#             canvas.print_png(buffer)
#             # 使用Image转换缓存数据
#             img = Image.open(buffer)
#             # 生成对应的array
#             img = np.asarray(img)
#             # 销毁
#             buffer.close()
#             # 写入gif
#             Writer.append_data(img)
```


text:
# Ordinary differential equations

## The Euler-Cromer scheme on a standard mesh

Consider the fundamental model problem for simple harmonic oscillations,

$$
\begin{equation}
u'' + \omega^2u = 0,\quad u(0)=I,\ u'(0)=0,
\label{stag:vib:ODE2nd} \tag{1}
\end{equation}
$$

where $\omega$ is the frequency of the oscillations (the exact solution
is $u(t)=I\cos\omega t$). This model can equivalently be
formulated as two first-order equations,

$$
\begin{equation}
v' = -\omega^2 u,
\label{stag:vib:EulerCromer:veq} \tag{2}
\end{equation}
$$

$$
\begin{equation}
u' = v\thinspace .
\label{stag:vib:EulerCromer:ueq} \tag{3}
\end{equation}
$$

The popular Euler-Cromer scheme for this $2\times 2$ system of ODEs
applies an explicit
forward difference in ([2](#stag:vib:EulerCromer:veq)) and a
backward difference in ([3](#stag:vib:EulerCromer:ueq)):

$$
\begin{equation}
\frac{v^{n+1}- v^n}{\Delta t} =- \omega^2u^{n},
\label{stag:vib:EulerCromer:dveq0} \tag{4}
\end{equation}
$$

$$
\begin{equation}
\frac{u^{n+1} - u^n}{\Delta t} = v^{n+1}\thinspace .
\label{stag:vib:EulerCromer:dueq0} \tag{5}
\end{equation}
$$

For a time domain $[0,T]$, we have introduced a mesh with points
$0=t_0 < t_1 < \cdots < t_n=T$. The most common case is a mesh with
uniform spacing $\Delta t$: $t_n=n\Delta t$. Then
$v^n$ is an approximation to $v(t)$ at mesh point $t_n$, and
$u^n$ is an approximation to $u(t)$ at the same point.
Note that the backward difference in ([7](#stag:vib:EulerCromer:dueq))
leads to an explicit updating formula for $u^{n+1}$
since $v^{n+1}$ is already computed:

$$
\begin{equation}
v^{n+1} = v^n -\Delta t \omega^2u^{n},
\label{stag:vib:EulerCromer:dveq} \tag{6}
\end{equation}
$$

$$
\begin{equation}
u^{n+1} = u^n + \Delta t v^{n+1}\thinspace .
\label{stag:vib:EulerCromer:dueq} \tag{7}
\end{equation}
$$

The Euler-Cromer scheme is equivalent with the standard second-order
accurate scheme for ([1](#stag:vib:ODE2nd)):

$$
\begin{equation}
u^{n+1} = 2u^n - u^{n-1} - \Delta t^2\omega^2 u^n,\ n=1,2,\ldots,
\label{stag:vib:dODE2nd} \tag{8}
\end{equation}
$$

but for the first time step, the method for ([1](#stag:vib:ODE2nd))
leads to

$$
\begin{equation}
u^1= u^0 - \frac{1}{2}\Delta t^2\omega^2u^0,
\label{_auto1} \tag{9}
\end{equation}
$$

while Euler-Cromer gives

$$
\begin{equation}
u^1= u^0 - \Delta t^2\omega^2u^0,
\label{_auto2} \tag{10}
\end{equation}
$$

which can be interpreted as a first-order, backward difference approximation
of $u'(0)=0$ combined with ([8](#stag:vib:dODE2nd)). At later time steps,
however, the alternating use of forward and backward differences in
([6](#stag:vib:EulerCromer:dveq))-([7](#stag:vib:EulerCromer:dueq))
leads to a method with error $\Oof{\Delta t^2}$.

# Exercises

## Exercise 1: Use the forward-backward scheme with quadratic damping

We consider the generalized model with quadratic damping, expressed
as a system of two first-order equations as in the section [vib:ode2:staggered](#vib:ode2:staggered):

$$
\begin{align*}
u^{\prime} &= v,\\
v' &= \frac{1}{m}\left( F(t) - \beta |v|v - s(u)\right)\thinspace .
\end{align*}
$$

However, contrary to what is done in the section [vib:ode2:staggered](#vib:ode2:staggered),
we want to apply the idea of a forward-backward discretization:
$u$ is marched forward by a one-sided Forward Euler scheme applied
to the first equation, and
thereafter $v$ can be marched forward by a Backward Euler scheme in the
second
equation.
Express the idea in operator notation and write out the
scheme. Unfortunately, the backward difference for the $v$ equation
creates a nonlinearity $|v^{n+1}|v^{n+1}$.  To linearize this
nonlinearity, use the known value $v^n$ inside the absolute value
factor, i.e., $|v^{n+1}|v^{n+1}\approx |v^n|v^{n+1}$.  Show that the
resulting scheme is equivalent to the one in the section [vib:ode2:staggered](#vib:ode2:staggered) for some time level $n\geq 1$.

What we learn from this exercise is that the first-order differences
and the linearization trick play together in "the right way" such that
the scheme is as good as when we (in the section [vib:ode2:staggered](#vib:ode2:staggered))
carefully apply centered differences and a geometric mean on a
staggered mesh to achieve second-order accuracy.
Filename: `vib_gen_bwdamping`.

# Partial differential equations

Examples on mesh functions on a staggered mesh in space.


text:
# Value at Risk and Conditional Value at Risk

We'll be looking at two important statistics of risk management when trading currencies: the Value at Risk (VaR) and the Conditional VaR (CVaR). As always, these notebooks are not concerned in teaching you the theoretical concepts in detail (there are enough webpages out there which will do a far better job). Instead, we brute force theory in your brain with raw mathematics and a hands-on practical example. Ahhh, the true data hacker way.

The main goal of this session is to show you why RISK is best seen as a flashing red light in your mind, deserving all your attention while trading. But what is risk? How do we estimate it? Between us, young padawan, we'll agree that risk is simply the amount of money you're about to lose. We need to somehow be able to estimate the risk we're about to take when we'll enter the market. Otherwise you'll go for coffee and come back with a brick in your window and your kidney gone. Stay safe out there kids, markets can move fast and banks are mafias.

The dataset used in this chapter is the EURUSD currency pair over a period of 4 years (2016 - 2019), sampled per minute. Because the file size is pretty large, I won't be sharing it on my slick hacker website. But if finding such data on the wild web is a problem, ask yourself: *'Who am I? am I a real worthy data hacker?'* Then poor an ice bucket over your head to reboot yourself, young one. If it were easy, everyone would do it!

## Long-tailed return distributions

Let's start off by looking at the amount a currency value can move in a given time: the return distribution.

```python
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.stats import norm, probplot

sns.set(style="ticks", context="talk")
plt.style.use("dark_background")

fn = "datasets/EURUSD_20162019_M1.csv"
data = pd.read_csv(fn, index_col = 0, parse_dates=True, dtype={'close':np.float64})#, 'ask_close':np.float64})
close = data.close * 10**4 #convert rate to Pips.
```

```python
#Set some plotting parameters based on the sample duration
for samplefreq in ['1min', '1H', '1D']:
    if samplefreq == '1min':
        close_sampled = close.copy()
        x_lim_range = 10
        bins = 800
        title_freq = 'minute'
    elif samplefreq == '1H':
        close_sampled = close.resample('1H').ohlc().close
        x_lim_range = 30
        bins = 500
        title_freq = 'hour'
    elif samplefreq == '1D':
        close_sampled = close.resample('1D').ohlc().close
        x_lim_range = 400
        bins = 100
        title_freq = 'day'

    #get returns
    close_return = close_sampled.diff()
    close_return = close_return.dropna(axis=0)

    #generate a normal distribution to compare
    mu = close_return.mean()
    sigma = close_return.std()
    x = np.linspace(-x_lim_range,x_lim_range,1000)
    y_norm = norm.pdf(x,loc=mu,scale=sigma)

    #PLOT histogram
    fs1 = 18
    fs2 = 16
    fs3 = 14

    fig, axes = plt.subplots(1,2,figsize=(16,6))
    sns.distplot(close_return, kde=True, bins=bins, label='data pdf', ax=axes[0])
    axes[0].set_xlim(-x_lim_range,x_lim_range)
    axes[0].set_title("EURUSD - 2016-2019 - Price movement per {}".format(title_freq), fontsize=fs1)
    axes[0].set_xlabel("Price movement [Pips]", fontsize=fs2)
    axes[0].set_ylabel("Probability", fontsize=fs2)
    #remove density line from seaborn
    axes[0].get_lines()[0].remove()
    #change fontsize of ticks
    axes[0].tick_params(axis='both', which='major', labelsize=fs3)
    #plot normal distributions
    axes[0].plot(x, y_norm, lw=2, alpha=0.6, label='normal pdf - same std as data', color='r')
    axes[0].legend(prop={'size': fs3})

    probplot(close_return, plot=axes[1])
    axes[1].set_title("QQ-plot", fontsize=fs1)
    axes[1].set_xlabel("Normally distributed theoretical quantiles", fontsize=fs2)
    axes[1].set_ylabel("Observed quantiles", fontsize=fs2)
    axes[1].tick_params(axis='both', which='major', labelsize=fs3)
    #plt.rcParams.update({'font.size': 16})

    plt.show()
```

B-e-a-utiful. If distributions don't turn you on than what does?
*The* most important information in these plots is that **intra-day price movements do not behave like a normal distribution.** Instead, they follow a **long-tailed distribution**. We can easily see that by looking at the first Quantile-Quantile (QQ) plot. To intuitively understand a QQ-plot, imagine taking the same amount of random samples from a normal distribution as we have in our dataset. If the dataset is normally distributed, the proportion of datasamples we have with small values in the middle compared to large values in the tails should be similar to the proportion of small to large values in an normal distribution. The per-minute data shows that way more samples are found in the tails of the distribution than we should expect had it been a normal distribution. If it were a normal distribution, the blue dots would have followed the red line (i.e. the quantiles would align).

This is veeeeery important, because many statistical approaches assume the noise in the market to be normally distributed! Starting with that assumption, statistical models instantly fail to provide a correct estimate of the risk! A model that assumes the price to follow a normal probability distribution will tell you there is about zero chance the price could drop 20 pips in one minute. However, the market shows that a 20 pip drop does happen, and it happens often!

What does all of this have to do with VaR and CVaR? Patience, young padawan. Eager you are to learn, I see. Very well. Stephen, bring in the math!

## Value-at-Risk (VaR)

Assume $\Delta V(t)$ equals the change in a currency pair, $t$ timesteps in the future, and $L(t)$ the corresponding lossfunction w.r.t. a position in the market. For a *long position*, $L(t) = -\Delta V(t)$ and for a *short position* $L(t) = \Delta V(t)$. The Value at Risk is then defined with corresponding probability $p$ as

\begin{equation}
    p = Pr[L(t) \geq VaR]
\end{equation}

In words: the probability that a position will result in a loss greater than or equal to the VaR over a timehorizon t is p. It follows that the VaR can be seen as the (1-p)-th quantile of the cumulatief distribution function of L(t):

\begin{equation}
     1-p = Pr[L(t) < VaR]
\end{equation}

In the context of the FX market, the VaR is mostly calculated for p = 0.01 or p = 0.05. The time horizon would then for instance be 1 hour or 1 day. The VaR corresponding to probability $1-p$ (last equation) can also be written as $VaR_{1-p}$.

## Conditional Value at Risk (CVaR)
**As the VaR only represents a quantile, it can lead to a significant underestimation of the potential loss if $\Delta V(t)$ follows a long-tailed distribution.** A better estimation of the potential loss is therefore the expected value of $L(t)$ given that the VaR is surpassed. This expected value is called the *Expected Shortfall* or *Conditional Value at Risk*:
\begin{equation}
    CVaR_{1-p} = \frac{1}{p} \int_{0}^{p} VaR_{1-q} dq
\end{equation}

## Demonstration

Let's have a look at the VaR and CVaR of the EURUSD pair sampled to t = 1 hour for p = 0.05. We'll assume we want to take a long position.

```python
#close_sampled = close.copy()
close_sampled = close.resample('1H').ohlc().close
close_return = close_sampled.diff()
close_return = close_return.dropna(axis=0)
p = 0.05

hist, edges = np.histogram(close_return, 1000)
pdf = hist / np.sum(hist)
cdf = np.cumsum(pdf)
var_index = np.where(cdf > p)[0][0]
var = edges[var_index]

cvar = np.sum(pdf[:var_index] * edges[:var_index] / np.sum(pdf[:var_index]))

print("*** EURUSD 2016-2019 return distribution / t=1hour / p=0.05 ***")
print("Value at Risk: {}\nConditional Value at Risk: {}".format(var,cvar))
```

    *** EURUSD 2016-2019 return distribution / t=1hour / p=0.05 ***
    Value at Risk: -15.28400000000218
    Conditional Value at Risk: -24.462376872966477

Intuitively: If we enter the market hoping it will go up and we come back in exactly one hour from now, there is a 5% chance the market will have dropped at least 15 pips. However, if the market drops, it may as well drop a lot further than 15 pips. GIVEN that the market drops at least 15 pips, the expected loss is 24 pips, a lot worse than the VaR.

Let's compare this result with a normal approximation of the market.

```python
#Check what Var and CVaR would be for a normal distribution with the same std.
mu = close_return.mean()
sigma = close_return.std()
var_norm = norm.ppf(0.05, loc=mu, scale=sigma) #ppf is the inverse cd
cvar_norm = norm(loc=mu, scale=sigma).expect(lambda x: x, lb=-1000, ub=var_norm, conditional=True)

print("*** Normal distribution (same mean and std as dataset) / p=0.05 ***")
print("Value at Risk: {}\nConditional Value at Risk: {}".format(var_norm,cvar_norm))
```

    *** Normal distribution (same mean and std as dataset) / p=0.05 ***
    Value at Risk: -17.186390585105393
    Conditional Value at Risk: -21.554713786264593

Because a normal distribution is fatter in the middle and dies out quicker in the tails, the inverse cdf of 0.05 is found at -17.186 pips while the true VaR is -15.284. The Value at Risk is therefore overestimated if the market were modelled with a normal distribution. The conditional expected value of the tail however, is only -21.554 while the true CVaR is larger (-24.462)! This effect of a long-tailed distribution only becomes worse if we had taken a smaller timehorizon.

Moral of the story:
1. Be careful with normal approximations
2. Risk is best quantified with the Conditional Value at Risk, rather than the VaR.

Note that we started out this chapter with a simple question: how do we estimate the risk we're exposing ourselves to? Some simple calculations later, we now have an actual number, the CVaR! It provides us with a quantified idea of the potential loss we should be able to take without selling our house and kids.

Alright, now it's time for the fun stuff: predicting the market! Onward!

max_stars_repo_name: serrHendrik/FX

max_issues_repo_name: serrHendrik/FX

max_forks_repo_name: serrHendrik/FX


text:
# 15.2. Solving equations and inequalities

https://ipython-books.github.io/152-solving-equations-and-inequalities/

## Ref

* SymPy Matrix guide at http://docs.sympy.org/latest/modules/matrices/matrices.html).
* Linear algebra on Wikipedia, at https://en.wikipedia.org/wiki/Linear_algebra#Further_reading
* Linear algebra on Wikibooks, at http://en.wikibooks.org/wiki/Linear_Algebra
* Linear algebra lectures on Awesome Math, at https://github.com/rossant/awesome-math/#linear-algebra

```python
from sympy import *
init_printing()
```

```python
var('x y z a')
```

```python
solve(x**2 - a, x)
```

```python
type(x)
```

    sympy.core.symbol.Symbol

```python
x = Symbol('x')
solve_univariate_inequality(x**2 > 4, x)
```

```python
x + 2*y + 1, x - 3*y - 2
```

```python
solve([x + 2*y + 1, x - 3*y - 2], x, y)
```

```python
solve([x + 2*y + 1, -x - 2*y - 1], x, y)
```

```python
solve([x**2 + y**2 - 1, x**2 - y**2 - S(1) / 2], x, y)
```

```python
var('a b c d u v')
```

```python
M = Matrix([[a, b, u], [c, d, v]])
M
```

$$\left[\begin{matrix}a & b & u\\c & d & v\end{matrix}\right]$$

```python
solve_linear_system(M, x, y)
```

```python
M[:2, :2]
```

$$\left[\begin{matrix}a & b\\c & d\end{matrix}\right]$$

```python
det(M[:2, :2])
```

```python

```


text:
# Scientific documents with $\LaTeX$ - another test

## Introduction

In your research, you will produce papers, reports and—very importantly—your thesis. These documents can be written using a WYSIWYG (What You See Is What You Get) editor (e.g., Word). However, an alternative especially suited for scientific publications is LaTeX. In LaTeX, the document is written  in a text file (`.tex`) with certain typesetting (tex) syntax. Text formatting is done using markups (like HTML). The file is then "compiled" (like source code of a programming language) into a file – typically in PDF.

### Why $\LaTeX$?

A number of reasons:

* The input is a small, portable text file
* LaTeX compilers are freely available for all OS'
* Exactly the same result on any computer (not true for Word, for example)
* LaTeX produces beautiful, professional looking docs
* Images are easy to embed and annotate
* Mathematical formulas (esp complex ones) are easy to write
* LaTeX is very stable – current version basically same since 1994! (9 major versions of MS Word since 1994 – with compatibility issues)
* LaTeX is free!
* You can focus on content, and not worry so much about formatting while writing
* An increasing number of Biology journals provide $\LaTeX$ templates, making formatting quicker.
* Referencing (bibliography) is easy (and can also be version controlled) and works with tools like Mendeley and Zotero
* Plenty of online support available – your question has probably already been answered
* You can integrate LaTeX into a workflow to auto-generate lengthy and complex documents (like your thesis).

---

**LaTeX documents scale up better then WYSIWYG editors.** Latex files for really large and complex documents (such as PhD Theses) are much easier to edit, manage and publish in specific formats (like pdf) than Word or Open Office documents.

---

### Limitations of $\LaTeX$

* It has a steeper learning curve.
* Can be difficult to manage revisions with multiple authors – especially if they don't use LaTeX! (Cue: Windows on a virtual machine!)
* Tracking changes are not available out of the box (but can be enabled using a suitable package)
* Typesetting tables can be a bit complex.
* Images and floats are easy to embed, and won't jump around like Word, but if you don't use the right package, they can be difficult to place where you want!

### Installing LaTeX

Type this in terminal:

```bash
sudo apt-get install texlive-full texlive-fonts-recommended texlive-pictures texlive-latex-extra imagemagick
```
It's a large installation, and will take some time.

We will use a text editor in this lecture, but you can use one of a number of dedicated editors (e.g., texmaker,
Gummi, TeXShop, etc.) There are also WYSIWYG frontends (e.g., Lyx, TeXmacs).

[Overleaf](https://www.overleaf.com/) is also very good (and works with git), especially for collaborating with non LaTeX-ers (your university may have a blanket license for the pro version).

## Key $\LaTeX$ features

### Environments

Environments are used to format blocks of text or graphics in a Latex document. They are delimited by an opening `\begin` and a closing `\end` tags (except for certain math environments). Everything inside will be formatted in a specific manner depending on the type of environment. For example, the code

```latex
\begin{center}
    Here is some text
\end{center}
```
Will produce "Here is some text" centered in the middle of the page.

The most commonly used Latex environments are:

|Environment| Purpose|
|:-|:-|
|`\begin{center} ... \end{center}`| Center the elements (works for text as well as graphics) |
|`\begin{itemize} ... \end{itemize}`| An itemized list  (default is bullet points)|
|`\begin{enumerate} ... \end{enumerate}`| An enumerated list (default is Arabic numerals) |
|`\begin{figure} ... \end{figure}`| For displaying a Figure |
|`\begin{table} ... \end{table}`| For displaying a table |
|`\begin{figure} ... \end{figure}`| For displaying a Figure |
|`\( ... \)`, `$ ... $`, or `\begin{math}...\end{math}`| For displaying an equation inline (as part of a sentence) |
|`\[ ... \]`, `$$ ... $$`, or `\begin{displaymath}...\end{displaymath}`| For producing an equation as a separate, display item (separate from the text)|
|`\begin{equation} ... \end{equation}`| For displaying a centered, numbered equation as a separate, display item|

In all of these environments, you can use modifier directives to the environment to tailor them. For example, in the `itemize` environment, the default is to create a list with bullets, but you can pick any symbol.

Below in your first Latex document example, you will see some examples of environments.

### Special characters

Some characters are "special" in Latex. These characters have a specific purpose, either inside a particular environment (e.g., table or equation), or both outside and inside an environment.

|Character| What it does|
|:-|:-|
|`#`| Used to reference arguments for a latex command; similar to the way `$` is an argument reference in [shell scripts](02-ShellScripting-Variables) |
|`$`| Used for opening or closing a mathematical equation or symbol; e.g.,`$y = mx + c$` gives $$y = mx + c$$|
|`%` | Comment character; everything from this symbol up to the end of line is ignored and will not appear in the final document |
|`&`| Alignment character; used to align columns in tables, and also equations in math environments|
|`_`| Subscript in math environments|
|`^`| Superscript in math environments|
|`{` and `}`| Use to group characters in math environments, and to enclose arguments in Latex commands|
|`~`| (equivalent to the command `\nobreakspace{}`) An "unbreakable" space, which can be used to add one or more "hard" spaces inside as well as outside math environments; for example, `$x  y$` gives $$x  y$$, but `$x~~y$` gives $$x~~y$$|
|`\`|  Indicates a LaTeX command, as in `\LaTeX` or `\maketitle` (or can be used to escape as special character - see below)|

#### Rendering special characters

If you want to actually reproduce these special characters in your document, you have to "escape" them by adding a *backslash* (`\`) in front of them. For example, writing `$\%$` produces the actual percentage symbol, $\%$.

### Latex commands

Every LaTeX command starts with a `\` . There are two types of commands:

* **Commands without arguments**: These commands are standalone, and do not take any additional arguments.
    * For example, in your first latex document above, the `\maketitle` command tells latex to render the title in the typeset document.
    * Another example: to render $\LaTeX$, you need the command `\LaTeX`

* **Commands with arguments**: These commands can (and often must) take arguments with curly brackets, which can be modified by including additional directives in square brackets before the main argument.
    * For example, in your first latex example above, the `\documentclass[12pt]{article}`
    * Another example: `\date{}` inserts the current date.

### Spaces and new lines

Note that:

* Several spaces in your text editor are treated as one space in the typeset document
* Several empty lines are treated as one empty line
* One empty line defines a new paragraph

### Typesetting math

There are two ways to display math

1. Inline (i.e., within the text).

2. Stand-alone, numbered equations and formulae.

For inline math, the "dollar" sign flanks the math to be typeset. For example, the code:

```
$\int_0^1 p^x (1-p)^y dp$
```

becomes $$\int_0^1 p^x (1-p)^y dp$$

For numbered equations, LaTeX provides the `equation` environment. For example,

```
\begin{equation}
    \int_0^1 \left(\ln \left( \frac{1}{x} \right)
    \right)^y dx = y!
\end{equation}
```

becomes

$$\int_0^1 \left(\ln \left( \frac{1}{x} \right) \right)^y dx = y!$$

### Document structure

Latex documents have a very specific structure in terms of the sequence in which certain elements must appear.

#### The start of the document

The first command is always `\documentclass[]{}` defining the type of document (e.g., `article, book, report, letter`).

Here, you can set several options. For example, to set size of text to 10 points and the letter paper size:

```latex
\documentclass[10pt,letterpaper]{article}
```

#### Defining packages

After having declared the type of document, you can specify special packages you want to use. Some particularly useful ones are:

|||
|:-|:-|
|`\usepackage{color}`| Use colors for text in your document|
|`\usepackage{amsmath,amssymb}`|Formats and commands for typesetting mathematical symbols and equations|
|`\usepackage{fancyhdr}`| Fine tune the formatting of headers and footers|
|`\usepackage{graphicx}`| Include figures in different formats: pdf, ps, eps, gif and jpeg|
|`\usepackage{listings}`| Typeset source code for different programming languages|
|`\usepackage{rotating}`| Allow rotation of tables and figures|
|`\usepackage{hyperref}`| Allow formatting of hyperlinks.
|`\usepackage{lineno}`| Allow line numbers|`_`| |
|

#### The main body

* Once you select the packages, you must start the main body of your document with `\begin{document}` and end it with `\end{document}`.

## A first LaTeX example

Let's try writing an example Latex document.

★ In your code editor, type the following and save it as a file called `FirstExample.tex` in a suitable location(e.g, in a `code` directory):

```latex

\documentclass[12pt]{article}

\title{A Simple Document}

\author{Your Name}

\date{}

\begin{document}
  \maketitle

  \begin{abstract}
    This paper analyzes a seminal equation in population biology.
  \end{abstract}

  \section{Introduction}
    Blah Blah

  \section{Materials \& Methods}

  A foundational equation of population biology is:

  \begin{equation}
    \frac{dN}{dt} = r N (1 - \frac{N}{K})
  \end{equation}

  It was first proposed by Verhulst in 1838 \cite{verhulst1838notice}.

  \bibliographystyle{plain}

  \bibliography{FirstBiblio}

\end{document}
```

```{note}
Look carefully at the way some of the elements such as special characters and environments are used in this first example document.
```

### Referencing and bibliography

Now, let's get a citation for the paper.

★ In the search box in Google Scholar type "verhulst population 1838"

The paper should be the only one (or the top one) to appear.

Click on the "Cite" icon (looks like two hollow commas) below the paper's title etc., and a small Cite window will appear. Click on "BibTeX" in the list of format options at the bottom, which should lead to a page with just the following text:

```bash
@article{verhulst1838notice,
  title={Notice sur la loi que la population suit dans son accroissement},
  author={Verhulst, Pierre-Fran{\c{c}}ois},
  journal={Corresp. Math. Phys.},
  volume={10},
  pages={113--126},
  year={1838}
}
```

Copy and paste this into a file called `FirstBiblio.bib`, saved in the same directory as `FirstExample.tex`

### Compiling the Latex document

Now we can create a `.pdf` file of the document.

★ In a terminal type (making sure you are the same directory where `FirstExample.tex` and `FirstBiblio.bib` are):

``` bash
 pdflatex FirstExample.tex
 bibtex FirstExample
 pdflatex FirstExample.tex
 pdflatex FirstExample.tex
```
This should produce the file `FirstExample.pdf`:

In the above bash script, we repeated the `pdflatex` command 3 times. Here's why:

* The first `pdflatex` run generates two files:`FirstExample.log` and `FirstExample.aux` (and an incomplete `.pdf`).
    * At this step, all cite{...} arguments info that bibtex needs are written into the `.aux` file.
* Then, the second `bibtex` command (followed by the filename without the `.tex` extension) results in bibtex reading the `.aux` file that was generated. It then produces two more files: `FirstExample.bbl` and `FirstExample.blg`
    * At this step, bibtex takes the citation info in the aux file and puts the relevant biblogrphic entries into the `.bbl` file (you can take a peek at all these files), formatted according to the instructions provided by the bibliography style that you have specified using `bibliographystyle{plain}`.
* The second `pdflatex` run updates `FirstExample.log` and `FirstExample.aux` (and a still-incomplete `.pdf` - the citations are not correctly formatted yet)
     * At this step, the reference list in the `.bbl` generated in the above step is included in the document, and the correct labels for the in-text `cite{...}` commands are written in `.aux` file (but the non in the actual pdf).
* The third and final `pdflatex` run then updates `FirstExample.log` and `FirstExample.aux` one last time, and now produces the complete `.pdf` file, with citations correctly formatted.
    * At this step, latex knows what the correct in-text citation labels are, and includes them in the pdf document.

Throughout all this, the `.log` file plays no role except to record info about how the commands are running.

PHEW! Why go through this repetitive sequence of commands? Well, "it is what it is" – $\LaTeX$, with all its advantages does have its quirks. The reason why it is this way, is probably that back then (Donald Knuth's PhD Thesis writing days – late 1950's to early 1960's), computers had *tiny* memories (RAMs), and writing files to disk and then reading them back in for the next step of the algorithm/program was the best (and only) way to go. Why has this not been fixed? I am not sure - keep an eye out, and it might well be (and then, raise an issue on TheMulQuaBio's [Github](https://github.com/mhasoba/TheMulQuaBio/issues)!)

Anyway, as such, you don't have to run these commands literally step by step, because you can create a bash script that does it for you, as we will now learn.

### A bash script to compile LaTeX

Let's write a useful little bash script to compile latex with bibtex.

★ Type the following script and call it `CompileLaTeX.sh` (you know where to put it!):

```bash
#!/bin/bash
pdflatex $1.tex
bibtex $1
pdflatex $1.tex
pdflatex $1.tex
evince $1.pdf &

## Cleanup
rm *.aux
rm *.log
rm *.bbl
rm *.blg
```
How do you run this script? The same as your previous bash scripts, so:

```latex
bash CompileLaTeX.sh FirstExample
```

### Exercise

Note that I have not written the `.tex` extension of `FirstExample` when feeding it to the in latex compilation bash script above. Make this bash script more convenient to use by allowing users to compile the script by using

```latex
bash CompileLaTeX.sh FirstExample.tex
```

## Some more $\LaTeX$ features and tips

Here are some more Latex features, and tips that might prove handy:

* LaTeX can render pretty much every mathematical symbol and operator that you can think of (plenty of lists and cheat-sheets online)
* Long documents can be split into separate `.tex` documents and combined using the `\input{}` command
* You can use bibliography managers such as [Mendeley](https://www.mendeley.com/?interaction_required=true) or [Zotero](https://www.zotero.org/) to export and maintain/update `.bib` files that are then ready to be used in a Latex document
* You can create new environments and commands, and create new ones in the preamble (which can also be kept as a  separate document and inserted using the `\input{}` command)

## Practicals

### First $\LaTeX$ example

Test `CompileLaTeX.sh` with `FirstExample.tex` and bring it under verson control under`week1` in your repository. Make sure that `CompileLaTeX.sh` will work if somebody else ran it from their computer using `FirstExample.tex` as an input.

### Practicals wrap-up

Make sure you have your `Week 1` directory organized with `Data`, `Sandbox` and `Code` with the necessary files and this week's (functional!) scripts in there. Every script should run without errors on my computer. This includes the five solutions (single-line commands you came up with) in `UnixPrac1.txt`.

*Commit and push every time you do some significant amount of coding work (after testing it), and then again before the given deadline (this will be announced in class).*

## Readings & Resources

### General

* [http://en.wikibooks.org/wiki/LaTeX/Introduction](http://en.wikibooks.org/wiki/LaTeX/Introduction)
* [The not so Short Introduction to LaTeX](https://ctan.org/tex-archive/info/lshort/english/)
* [The Visual LaTeX FAQ: sometimes it is difficult to describe what you want to do!](http://mirror.las.iastate.edu/tex-archive/info/visualFAQ/visualFAQ.pdf)
* [The Overleaf knowledge base](https://www.overleaf.com/learn), including
    * [Learn LaTeX in 30 minutes](https://www.overleaf.com/learn/latex/Learn_LaTeX_in_30_minutes)
    * [Presentations in LaTeX](https://www.overleaf.com/learn/latex/Beamer_Presentations:_A_Tutorial_for_Beginners_(Part_1)—Getting_Started)
    * [Bibliographies in LaTeX](https://www.overleaf.com/learn/latex/Bibliography_management_with_bibtex)

### LaTeX Templates
* There are lots of LaTeX templates online, such for typesetting theses from particular institutions, or papers for a specific journal. There are some examples the `TheMulQuaBio` repo (under `code`).
* The [Overleaf templates](https://www.overleaf.com/latex/templates) are extensive
    * These includes many [Imperial College Dissertation templates](https://www.overleaf.com/latex/templates?addsearch=imperial%20college)).

### $\LaTeX$ Tables
* [$\LaTeX$ table generator](http://www.tablesgenerator.com/)

max_stars_repo_name: mhasoba/TheMulQuaBio

max_issues_repo_name: mhasoba/TheMulQuaBio

max_forks_repo_name: mhasoba/TheMulQuaBio


text:
# Problema de selección de portafolio con preferencias media-varianza

En la clase pasada hablamos acerca de:
- preferencias,
- funciones de utilidad,
- la actitud de los inversionistas de cara al riesgo,
- la aversión al riesgo, entre otros.

Todas ellas son piezas que necesitamos para responder la pregunta de ¿cómo un inversionista toma la decisión óptima de selección de portafolio?

En esta clase al fin estamos listos para ensamblar estos conceptos y escribir el problema de selección de portafolios.

En el camino aprenderemos acerca del concepto de **utilidad esperada**, que nos permite trabajar con incertidumbre en el modelado económico (una de las ideas más importantes en economía). Esta idea tiene más de 60 años, y básicamente dice que los individuos, cuando están de cara a incertidumbre, maximizan el valor esperado de su utilidad (solo cierto si somos homo economicus).

Además del concepto de utilidad esperada, aprenderemos acerca de **preferencias media-varianza**. Es decir, supondremos que los inversionistas toman decisiones basados en un tipo particular de preferencias.

Con lo anterior, estableceremos el problema de selección de portafolios.

**Objetivos:**

- ¿Qué es utilidad esperada?
- ¿Qué son las preferencias media-varianza?
- Funciones de utilidad media-varianza.
- Enunciado y solución del problema básico de selección de portafolio.

*Referencia:*
- Notas del curso "Portfolio Selection and Risk Management", Rice University, disponible en Coursera.
___

## 1. Utilidad esperada
- Básicamente combina las probabilidades de los resultados con cómo los inversionistas se sienten con dichos resultados.
- En otras palabras, la utilidad esperada multiplica la probabilidad de suceso de un evento con la utilidad que genera dicho evento.

Recordemos que las *funciones de utilidad* permiten a los inversionistas expresar cómo se sienten con los resultados, especialmente en los malos ratos.

Entonces la *utilidad esperada* es una herramienta que nos permite cuantificar cómo nos sentimos en nuestros malos momentos económicos, capturando el riesgo con la probabilidad de ocurrencia de dichos malos momentos.

Dado este marco de trabajo, cualquier decisión se puede escribir como la maximización de la utilidad esperada:
\begin{align}
\max_{\theta} & \quad E[U(W)], \\
\end{align}
mediante la escogencia de cierta variable $\theta$ (gastos, planes de ahorro, compra de activos, planes de producción, etc.).

Para nuestros propósitos, la variable de decisión serán los pesos o ponderaciones del portafolio.

Adicionalmente, en el contexto de la decisión de distribución de la riqueza entre activos, el problema de maximización tendrá comúnmente las siguientes restricciones:
- universo de inversión,
- posición en los activos dados.

**Ejemplo.**

Supongamos que un inversionista debe determinar la composición óptima de su portafolio, que contiene activos y bonos. Supongamos que son los únicos instrumentos disponibles.

Sean:
- $w_s$: peso o ponderación de activos en el portafolio,
- $w_b$: peso o ponderación de bonos en el portafolio,
- $r_s$: rendimiento de los activos, y
- $r_b$: rendimiento de los bonos.

De manera que podemos escribir el problema de selección de portafolios como la maximización de la utilidad esperade de nuestra riqueza futura, la cual dependerá de nuestros rendimientos:

\begin{align}
\max_{w_s,w_b} &\quad E[U(W)]\\
\text{s. a.}   &\quad W=W_0(1+w_sr_s+w_br_b)\\
               &\quad w_s+w_b=1\\
               &\quad w_s,w_b\geq0
\end{align}

Preguntas:
- ¿Qué significan las restricciones?
- Ya que tenemos planteado este problema básico, ¿qué haría falta para empezar a resolverlo?
___

## 2. Preferencias media-varianza

### 2.1. Utilidad media-varianza

Entonces, ¿qué funciones de utilidad deberíamos de usar en este problema de selección de portafolios?

- La respuesta es: **preferencias media-varianza**.
- Éstas serán representadas en términos de funciones de utilidad como: **utilidad media-varianza**.

Usamos la *utilidad media-varianza* en el problema de selección de portafolios dado que ésta decribe el "trade-off" entre riesgo y rendimiento que enfrentan los inversionistas. La *utilidad media-varianza* está dada por la siguiente expresión:

$$U=E[r_p]-\frac{1}{2}\gamma\sigma_p^2,$$

donde
- $E[r_p]$ es el rendimiento esperado del portafolio,
- $\sigma_p^2$ es la varianza del portafolio, y
- $\gamma$ es el coeficiente de aversión al riesgo.

#### Intuición acerca de la función de utilidad media-varianza:
- Sólo se preocupa por medias :) y varianzas :(.
- Incrementa con: rendimiento esperado del portafolio.
- Decrece con: varianza del portafolio.
- Malos tiempos: rendimientos son bajos y las volatilidades son altas.
- Conecta bastante bien con la teoría moderna de portafolios, la cual caracteriza los rendimientos con medias y varianzas únicamente.
- Criticada por su limitación: supone que los inversionistas sólo se preocupan por medias y varianzas.

### 2.2. Curvas de indiferencia

*¿Recuerdan las curvas de nivel que se ven en cálculo de varias variables?*
- Bien, acá nos servirán para representar la utilidad media-varianza gráficamente.
- En el contexto de utilidad media-varianza, las curvas de nivel se llaman **curvas de indiferencia**.

Dados ciertos niveles de utilidad $U_1>U_2>U_3$, las curvas de indiferencia relativas a estos niveles de utilidad, son los lugares geométricos en el espacio de rendimiento esperado vs. volatilidad representados por las siguientes expresiones

$$U_1=E[r_p]-\frac{1}{2}\gamma\sigma_p^2\Rightarrow E[r_p]=\frac{1}{2}\gamma\sigma_p^2+U_1,$$

$$U_2=E[r_p]-\frac{1}{2}\gamma\sigma_p^2\Rightarrow E[r_p]=\frac{1}{2}\gamma\sigma_p^2+U_2,$$

$$U_3=E[r_p]-\frac{1}{2}\gamma\sigma_p^2\Rightarrow E[r_p]=\frac{1}{2}\gamma\sigma_p^2+U_3.$$

**Gráficamente**

```python
# Importar numpy y pyplot
import numpy as np
from matplotlib import pyplot as plt
%matplotlib inline
```

```python
# Coeficiente de aversión al riesgo (entre 1 y 10 comúnmente)
g = 7
# Niveles de utilidad
U1, U2, U3 = 0.35, 0.2, 0.05
# Vector de volatilidades (sugerido 0%-60%)
sp = np.linspace(0, 0.6)
# Curvas de indiferencia
Erp1 = 1 / 2 * g * sp**2 + U1
Erp2 = 1 / 2 * g * sp**2 + U2
Erp3 = 1 / 2 * g * sp**2 + U3
```

```python
# Gráfica
plt.figure(figsize=(6, 4))
plt.plot(sp, Erp1, label=f'Curva de indiferencia $U_1$={U1}')
plt.plot(sp, Erp2, label=f'Curva de indiferencia $U_2$={U2}')
plt.plot(sp, Erp3, label=f'Curva de indiferencia $U_3$={U3}')
plt.xlabel('Volatilidad $\sigma$')
plt.ylabel('Rendimiento esperado $E[r]$')
plt.grid()
plt.legend(loc='upper left', bbox_to_anchor=(1.05, 1))
```

Bueno, ¿y porqué se llaman curvas de indiferencia?, ¿qué representa una cuva de indiferencia?

- Porque sobre una misma curva el nivel de utilidad es el mismo (es indiferente).
- Son todas las combinaciones de riesgo y rendimiento que producen un mismo nivel de utilidad.

Volviendo al problema de selección de portafolios, queremos la utilidad más alta.
- ¿Cuál de las anteriores curvas de indiferencia corresponde a la utilidad más alta?
- Intuitivamente, ¿porqué?
- Curvas de indiferencia para niveles de utilidad más altos, estarán...

Notamos además que las anteriores curvas de indiferencia son *paralelas* una con otra. Claro, las dibujamos con el mismo coeficiente de aversión al riesgo.

$$\frac{d E[r_p]}{d\sigma_p}=\frac{d }{d\sigma_p}\left[\frac{1}{2}\gamma\sigma_p^2+U\right]=\gamma\sigma_p.$$

¿Cómo cambian estas curvas para coeficientes de aversión al riesgo más altos?

```python
# Coeficientes de aversión al riesgo (entre 1 y 10 comúnmente)
g1, g2, g3 = 3, 5, 7
# Nivel de utilidad
U = 0.2
# Vector de volatilidades (sugerido 0%-60%)
sp = np.linspace(0, 0.6)
# Curvas de indiferencia
Erp1 = 1 / 2 * g1 * sp**2 + U
Erp2 = 1 / 2 * g2 * sp**2 + U
Erp3 = 1 / 2 * g3 * sp**2 + U
```

```python
# Gráfica
plt.figure(figsize=(6, 4))
plt.plot(sp, Erp1, label=f'Curva de indiferencia $\gamma_1$={g1}')
plt.plot(sp, Erp2, label=f'Curva de indiferencia $\gamma_2$={g2}')
plt.plot(sp, Erp3, label=f'Curva de indiferencia $\gamma_3$={g3}')
plt.xlabel('Volatilidad $\sigma$')
plt.ylabel('Rendimiento esperado $E[r]$')
plt.grid()
plt.legend(loc='upper left', bbox_to_anchor=(1.05, 1))
```

¿Cómo interpretamos las anteriores gráficas?, ¿qué pasa con las personas más aversas al riesgo?
- Se puede ver de dos maneras: para un mismo nivel de rendimiento esperado, una persona más aversa al riesgo soporta un nivel menor de riesgo; equivalentemente, para un mismo nivel de riesgo, una persona más aversa al riesgo requerirá un nivel de rendimiento esperado más alto.

Con todo lo anterior, el problema de selección de portafolios se puede plantear como *encontrar la curva de indeferencia más alta dado el conjunto de oportunidades de inversión y restricciones*.

## 3. Problema de selección de portafolios: una ilustración

Ahora ilustraremos el problema de selección de portafolios con algunos datos.
- Por ahora solo queremos ilustrar gráficamente cómo se resuelve este problema. Trabajar en la intuición.
- En las siguientes dos clases nos enfocaremos en cómo resolverlo analíticamente.

Acá tenemos el rendimiento medio anual y la volatilidad para dos instrumentos usando datos de EU: instrumentos de deuda (bonos) y acciones. Supondremos que el inversionista solo puede invertir en estas dos clases de instrumentos.

```python
# Importamos pandas
import pandas as pd
```

```python
# Datos
data = pd.DataFrame(index=['Stocks','Bonds', 'CorrSB'], columns=['Mean', 'Std'])
data['Mean'] = [0.119, 0.0591, 0.113]
data['Std'] = [0.1915, 0.0833, None]
data
```

      Mean
      Std

      Stocks
      0.1190
      0.1915

      Bonds
      0.0591
      0.0833

      CorrSB
      0.1130
      NaN

Entonces, ¿cuál es la distribución de riqueza óptima?, o más bien, ¿cuál es la composición óptima del portafolio para un inversionista dado su nivel de aversión al riesgo?

**Primero.** Recordamos que, para dos activos, podemos trazar la frontera de mínima varianza tomando todas las posibles combinaciones de los dos activos.

De nuevo, sean:
- $w_s=w$: peso o ponderación de activos en el portafolio,
- $w_b=1-w$: peso o ponderación de bonos en el portafolio,
- $r_s$: rendimiento de los activos, y
- $r_b$: rendimiento de los bonos.

Entonces

$$E[r_p]=wE[r_{s}]+(1-w)E[r_b]$$

$$\sigma_p^2=w^2\sigma_{s}^2+(1-w)^2\sigma_b^2+2w(1-w)\rho_{s,b}\sigma_s\sigma_b$$

```python
# Vector de w variando entre 0 y 1 con n pasos
n = 101
w = np.linspace(0, 1, n)
# Rendimientos esperados individuales
Ers = data.loc['Stocks', 'Mean']
Erb = data.loc['Bonds', 'Mean']
# Volatilidades individuales
ss = data.loc['Stocks', 'Std']
sb = data.loc['Bonds', 'Std']
# Correlacion
rsb = data.loc['CorrSB', 'Mean']
```

```python
# Crear un DataFrame cuyas columnas sean rendimiento
# y volatilidad del portafolio para cada una de las w
# generadas
portafolios = pd.DataFrame(data={'w': w,
                                 'Media': w * Ers + (1 - w) * Erb,
                                 'Vol': (w**2 * ss**2 + (1 - w)**2 * sb**2 + 2 * w * (1 - w) * ss * sb * rsb)**0.5
                                })
portafolios
```

      w
      Media
      Vol

      0
      0.00
      0.059100
      0.083300

      1
      0.01
      0.059699
      0.082705

      2
      0.02
      0.060298
      0.082155

      3
      0.03
      0.060897
      0.081650

      4
      0.04
      0.061496
      0.081191

      ...
      ...
      ...
      ...

      96
      0.96
      0.116604
      0.184246

      97
      0.97
      0.117203
      0.186054

      98
      0.98
      0.117802
      0.187866

      99
      0.99
      0.118401
      0.189681

      100
      1.00
      0.119000
      0.191500

101 rows × 3 columns

```python
# Gráfica
plt.figure(figsize=(6, 4))
plt.plot(portafolios['Vol'], portafolios['Media'], label='Portafolios')
plt.plot(sb, Erb, 'or', label='Bonds')
plt.plot(ss, Ers, 'ob', label='Stocks')
plt.xlabel('Volatilidad $\sigma$')
plt.ylabel('Rendimiento esperado $E[r]$')
plt.grid()
plt.legend(loc='upper left', bbox_to_anchor=(1.05, 1))
```

**Segundo.** Graficamos en la misma ventana, curvas de indiferencia.

```python
# Niveles de utilidad
U1, U2, U3 = 0.05, 0.06, 0.0723
# Coeficiente de aversión al riesgo
g = 3
# Curvas de indiferencia
sp = np.linspace(0.06, 0.2)
Erp1 = 1 / 2 * g * sp**2 + U1
Erp2 = 1 / 2 * g * sp**2 + U2
Erp3 = 1 / 2 * g * sp**2 + U3
```

```python
# Gráfica
plt.figure(figsize=(6, 4))
plt.plot(portafolios['Vol'], portafolios['Media'], label='Portafolios')
plt.plot(sb, Erb, 'or', label='Bonds')
plt.plot(ss, Ers, 'ob', label='Stocks')
plt.plot(sp, Erp1, label=f'Curva de indiferencia $U_1$={U1}')
plt.plot(sp, Erp2, label=f'Curva de indiferencia $U_2$={U2}')
plt.plot(sp, Erp3, label=f'Curva de indiferencia $U_3$={U3}')
plt.xlabel('Volatilidad $\sigma$')
plt.ylabel('Rendimiento esperado $E[r]$')
plt.grid()
plt.legend(loc='upper left', bbox_to_anchor=(1.05, 1))
#plt.axis([0.09, 0.11, 0.08, 0.09])
```

**Tercero.** La elección óptima está dada por la curva de indiferencia para el nivel de utilidad más alto que es tangente a la frontera media-varianza.
- Claramente, esta selección depende del coeficiente de aversión al riesgo.

```python
# Gráfica con zoom
# Gráfica
plt.figure(figsize=(6, 4))
plt.plot(portafolios['Vol'], portafolios['Media'], label='Portafolios')
plt.plot(sb, Erb, 'or', label='Bonds')
plt.plot(ss, Ers, 'ob', label='Stocks')
plt.plot(sp, Erp1, label=f'Curva de indiferencia $U_1$={U1}')
plt.plot(sp, Erp2, label=f'Curva de indiferencia $U_2$={U2}')
plt.plot(sp, Erp3, label=f'Curva de indiferencia $U_3$={U3}')
plt.xlabel('Volatilidad $\sigma$')
plt.ylabel('Rendimiento esperado $E[r]$')
plt.grid()
plt.legend(loc='upper left', bbox_to_anchor=(1.05, 1))
plt.axis([0.12, 0.14, 0.09, 0.105])
```

# Anuncios parroquiales

## 1. Quiz la siguiente clase.

Created with Jupyter by Esteban Jiménez Rodríguez.


text:
# Full-Waveform Inversion (FWI)

This notebook highlights various aspects of seismic inversion based on Devito operators. In this example we aim to illustrate the core ideas behind seismic inversion, where we create an image of the subsurface from field recorded data. This tutorial follows on `03_propagators-acoustic.ipynb` notebook.

## Inversion requirement

Seismic inversion relies on two known parameters:

- **Field data** - or also called **recorded data**. This is a shot record corresponding to the true velocity model. In practice this data is acquired as described in the first tutorial. In order to simplify this tutorial we will fake field data by modelling it with the true velocity model.

- **Initial velocity model**. This is a velocity model that has been obtained by processing the field data. This model is a rough and very smooth estimate of the velocity as an initial estimate for the inversion. This is a necessary requirement for any optimization (method).

## Inversion computational setup

We will introduce the gradient operator. This operator corresponds to the imaging condition introduced in the previous tutorial with some minor modifications that are defined by the objective function (also referred to in the tutorial series as the *functional*, *f*) and its gradient, *g*. We will also define these two terms in the tutorial.

## Notes on the operators

In `03_propagators-acoustic.ipynb` we described in detail the creation of a forward modelling operator -- here we will call a wrapper function to such a model instead. This wrapper contains all necessary operators for seismic modeling, imaging and inversion. Operators introduced for the first time in this tutorial will be properly described.

```python
import numpy as np
%matplotlib inline

from devito import configuration
configuration['log-level'] = 'WARNING'
```

## Computational considerations

As we will see, FWI is extremely computationally demanding (more so than, e.g., RTM). To keep this tutorial as lightwight as possible we therefore again use a very small demonstration model. We also define here a few parameters for the final example runs that can be changed to modify the overall runtime of the tutorial.

```python
nshots = 9  # Number of shots to create gradient from
nreceivers = 101  # Number of receiver locations per shot
fwi_iterations = 5  # Number of outer FWI iterations
```

# True and smooth velocity models

We will use a very simple model domain, consisting of a circle within a 2D domain. We will again use the "true" model to generate our synthetic shot data and use a "smooth" model as our initial guess. In this case the smooth model is very smooth indeed - it is simply a constant background velocity without any features.

```python
#NBVAL_IGNORE_OUTPUT
from examples.seismic import demo_model, plot_velocity, plot_perturbation

# Define true and initial model
shape = (101, 101)  # Number of grid point (nx, nz)
spacing = (10., 10.)  # Grid spacing in m. The domain size is now 1km by 1km
origin = (0., 0.)  # Need origin to define relative source and receiver locations

model = demo_model('circle-isotropic', vp_circle=3.0, vp_background=2.5,
                    origin=origin, shape=shape, spacing=spacing, nbl=40)

model0 = demo_model('circle-isotropic', vp_circle=2.5, vp_background=2.5,
                     origin=origin, shape=shape, spacing=spacing, nbl=40,
                     grid = model.grid)

plot_velocity(model)
plot_velocity(model0)
plot_perturbation(model0, model)
```

```python
assert model.grid == model0.grid
assert model.vp.grid == model0.vp.grid
```

## Acquisition geometry

In this tutorial, we will use the easiest case for inversion, namely a transmission experiment. The sources are located on one side of the model and the receivers on the other side. This allows to record most of the information necessary for inversion, as reflections usually lead to poor inversion results.

```python
#NBVAL_IGNORE_OUTPUT
# Define acquisition geometry: source
from examples.seismic import AcquisitionGeometry

t0 = 0.
tn = 1000.
f0 = 0.010
# First, position source centrally in all dimensions, then set depth
src_coordinates = np.empty((1, 2))
src_coordinates[0, :] = np.array(model.domain_size) * .5
src_coordinates[0, 0] = 20.  # Depth is 20m

# Define acquisition geometry: receivers

# Initialize receivers for synthetic and imaging data
rec_coordinates = np.empty((nreceivers, 2))
rec_coordinates[:, 1] = np.linspace(0, model.domain_size[0], num=nreceivers)
rec_coordinates[:, 0] = 980.

# Geometry

geometry = AcquisitionGeometry(model, rec_coordinates, src_coordinates, t0, tn, f0=f0, src_type='Ricker')
# We can plot the time signature to see the wavelet
geometry.src.show()
```

```python
# Plot acquisition geometry
plot_velocity(model, source=geometry.src_positions,
              receiver=geometry.rec_positions[::4, :])
```

## True and smooth data

We can generate shot records for the true and smoothed initial velocity models, since the difference between them will again form the basis of our imaging procedure.

```python
# Compute synthetic data with forward operator
from examples.seismic.acoustic import AcousticWaveSolver

solver = AcousticWaveSolver(model, geometry, space_order=4)
true_d, _, _ = solver.forward(vp=model.vp)
```

```python
# Compute initial data with forward operator
smooth_d, _, _ = solver.forward(vp=model0.vp)
```

```python
#NBVAL_IGNORE_OUTPUT
from examples.seismic import plot_shotrecord

# Plot shot record for true and smooth velocity model and the difference
plot_shotrecord(true_d.data, model, t0, tn)
plot_shotrecord(smooth_d.data, model, t0, tn)
plot_shotrecord(smooth_d.data - true_d.data, model, t0, tn)
```

# Full-Waveform Inversion

## Formulation

Full-waveform inversion (FWI) aims to invert an accurate model of the discrete wave velocity, $\mathbf{c}$, or equivalently the square slowness of the wave, $\mathbf{m} = \frac{1}{\mathbf{c}^2}$, from a given set of measurements of the pressure wavefield $\mathbf{u}$. This can be expressed as the following optimization problem [1, 2]:

\begin{aligned}
	\mathop{\hbox{minimize}}_{\mathbf{m}} \Phi_s(\mathbf{m})&=\frac{1}{2}\left\lVert\mathbf{P}_r
	\mathbf{u} - \mathbf{d}\right\rVert_2^2 \\
	\mathbf{u} &= \mathbf{A}(\mathbf{m})^{-1} \mathbf{P}_s^T \mathbf{q}_s,
\end{aligned}

where $\mathbf{P}_r$ is the sampling operator at the receiver locations, $\mathbf{P}_s^T$ is the injection operator at the source locations, $\mathbf{A}(\mathbf{m})$ is the operator representing the discretized wave equation matrix, $\mathbf{u}$ is the discrete synthetic pressure wavefield, $\mathbf{q}_s$ is the corresponding pressure source and $\mathbf{d}$ is the measured data. It is worth noting that $\mathbf{m}$ is the unknown in this formulation and that multiple implementations of the wave equation operator $\mathbf{A}(\mathbf{m})$ are possible.

We have already defined a concrete solver scheme for $\mathbf{A}(\mathbf{m})$ in the first tutorial, including appropriate implementations of the sampling operator $\mathbf{P}_r$ and source term $\mathbf{q}_s$.

To solve this optimization problem using a gradient-based method, we use the
adjoint-state method to evaluate the gradient $\nabla\Phi_s(\mathbf{m})$:

\begin{align}
 \nabla\Phi_s(\mathbf{m})=\sum_{\mathbf{t} =1}^{n_t}\mathbf{u}[\mathbf{t}] \mathbf{v}_{tt}[\mathbf{t}] =\mathbf{J}^T\delta\mathbf{d}_s,
\end{align}

where $n_t$ is the number of computational time steps, $\delta\mathbf{d}_s = \left(\mathbf{P}_r \mathbf{u} - \mathbf{d} \right)$ is the data residual (difference between the measured data and the modelled data), $\mathbf{J}$ is the Jacobian operator and $\mathbf{v}_{tt}$ is the second-order time derivative of the adjoint wavefield solving:

\begin{align}
 \mathbf{A}^T(\mathbf{m}) \mathbf{v} = \mathbf{P}_r^T \delta\mathbf{d}.
\end{align}

We see that the gradient of the FWI function is the previously defined imaging condition with an extra second-order time derivative. We will therefore reuse the operators defined previously inside a Devito wrapper.

## FWI gradient operator

To compute a single gradient $\nabla\Phi_s(\mathbf{m})$ in our optimization workflow we again use `solver.forward` to compute the entire forward wavefield $\mathbf{u}$ and a similar pre-defined gradient operator to compute the adjoint wavefield `v`. The gradient operator provided by our `solver` utility also computes the correlation between the wavefields, allowing us to encode a similar procedure to the previous imaging tutorial as our gradient calculation:

- Simulate the forward wavefield with the background velocity model to get the synthetic data and save the full wavefield $\mathbf{u}$
- Compute the data residual
- Back-propagate the data residual and compute on the fly the gradient contribution at each time step.

This procedure is applied to multiple source positions and summed to obtain a gradient image of the subsurface. We again prepare the source locations for each shot and visualize them, before defining a single gradient computation over a number of shots as a single function.

```python
#NBVAL_IGNORE_OUTPUT

# Prepare the varying source locations sources
source_locations = np.empty((nshots, 2), dtype=np.float32)
source_locations[:, 0] = 30.
source_locations[:, 1] = np.linspace(0., 1000, num=nshots)

plot_velocity(model, source=source_locations)
```

```python
from devito import Eq, Operator

# Computes the residual between observed and synthetic data into the residual
def compute_residual(residual, dobs, dsyn):
    if residual.grid.distributor.is_parallel:
        # If we run with MPI, we have to compute the residual via an operator
        # First make sure we can take the difference and that receivers are at the
        # same position
        assert np.allclose(dobs.coordinates.data[:], dsyn.coordinates.data)
        assert np.allclose(residual.coordinates.data[:], dsyn.coordinates.data)
        # Create a difference operator
        diff_eq = Eq(residual, dsyn.subs({dsyn.dimensions[-1]: residual.dimensions[-1]}) -
                               dobs.subs({dobs.dimensions[-1]: residual.dimensions[-1]}))
        Operator(diff_eq)()
    else:
        # A simple data difference is enough in serial
        residual.data[:] = dsyn.data[:] - dobs.data[:]

    return residual
```

```python
# Create FWI gradient kernel
from devito import Function, TimeFunction, norm
from examples.seismic import Receiver

import scipy
def fwi_gradient(vp_in):
    # Create symbols to hold the gradient
    grad = Function(name="grad", grid=model.grid)
    # Create placeholders for the data residual and data
    residual = Receiver(name='residual', grid=model.grid,
                        time_range=geometry.time_axis,
                        coordinates=geometry.rec_positions)
    d_obs = Receiver(name='d_obs', grid=model.grid,
                     time_range=geometry.time_axis,
                     coordinates=geometry.rec_positions)
    d_syn = Receiver(name='d_syn', grid=model.grid,
                     time_range=geometry.time_axis,
                     coordinates=geometry.rec_positions)
    objective = 0.
    for i in range(nshots):
        # Update source location
        geometry.src_positions[0, :] = source_locations[i, :]

        # Generate synthetic data from true model
        _, _, _ = solver.forward(vp=model.vp, rec=d_obs)

        # Compute smooth data and full forward wavefield u0
        _, u0, _ = solver.forward(vp=vp_in, save=True, rec=d_syn)

        # Compute gradient from data residual and update objective function
        compute_residual(residual, d_obs, d_syn)

        objective += .5*norm(residual)**2
        solver.gradient(rec=residual, u=u0, vp=vp_in, grad=grad)

    return objective, grad
```

Having defined our FWI gradient procedure we can compute the initial iteration from our starting model. This allows us to visualize the gradient alongside the model perturbation and the effect of the gradient update on the model.

```python
# Compute gradient of initial model
ff, update = fwi_gradient(model0.vp)
assert np.isclose(ff, 57283, rtol=1e0)
```

```python
#NBVAL_IGNORE_OUTPUT
from devito import mmax
from examples.seismic import plot_image

# Plot the FWI gradient
plot_image(-update.data, vmin=-1e4, vmax=1e4, cmap="jet")

# Plot the difference between the true and initial model.
# This is not known in practice as only the initial model is provided.
plot_image(model0.vp.data - model.vp.data, vmin=-1e-1, vmax=1e-1, cmap="jet")

# Show what the update does to the model
alpha = .5 / mmax(update)
plot_image(model0.vp.data + alpha*update.data, vmin=2.5, vmax=3.0, cmap="jet")
```

We see that the gradient and the true perturbation have the same sign, therefore, with an appropriate scaling factor, we will update the model in the correct direction.

```python
from sympy import Min, Max
# Define bounding box constraints on the solution.
def update_with_box(vp, alpha, dm, vmin=2.0, vmax=3.5):
    """
    Apply gradient update in-place to vp with box constraint

    Notes:
    ------
    For more advanced algorithm, one will need to gather the non-distributed
    velocity array to apply constrains and such.
    """
    update = vp + alpha * dm
    update_eq = Eq(vp, Max(Min(update, vmax), vmin))
    Operator(update_eq)()
```

```python
#NBVAL_SKIP

from devito import mmax

# Run FWI with gradient descent
history = np.zeros((fwi_iterations, 1))
for i in range(0, fwi_iterations):
    # Compute the functional value and gradient for the current
    # model estimate
    phi, direction = fwi_gradient(model0.vp)

    # Store the history of the functional values
    history[i] = phi

    # Artificial Step length for gradient descent
    # In practice this would be replaced by a Linesearch (Wolfe, ...)
    # that would guarantee functional decrease Phi(m-alpha g) <= epsilon Phi(m)
    # where epsilon is a minimum decrease constant
    alpha = .05 / mmax(direction)

    # Update the model estimate and enforce minimum/maximum values
    update_with_box(model0.vp , alpha , direction)

    # Log the progress made
    print('Objective value is %f at iteration %d' % (phi, i+1))
```

    Objective value is 39291.434160 at iteration 1
    Objective value is 24504.771842 at iteration 2
    Objective value is 14385.612613 at iteration 3
    Objective value is 7907.016246 at iteration 4
    Objective value is 3959.785426 at iteration 5

```python
#NBVAL_IGNORE_OUTPUT

# Plot inverted velocity model
plot_velocity(model0)
```

```python
#NBVAL_SKIP
import matplotlib.pyplot as plt

# Plot objective function decrease
plt.figure()
plt.loglog(history)
plt.xlabel('Iteration number')
plt.ylabel('Misift value Phi')
plt.title('Convergence')
plt.show()
```

## References

[1] _Virieux, J. and Operto, S.: An overview of full-waveform inversion in exploration geophysics, GEOPHYSICS, 74, WCC1–WCC26, doi:10.1190/1.3238367, http://library.seg.org/doi/abs/10.1190/1.3238367, 2009._

[2] _Haber, E., Chung, M., and Herrmann, F. J.: An effective method for parameter estimation with PDE constraints with multiple right hand sides, SIAM Journal on Optimization, 22, http://dx.doi.org/10.1137/11081126X, 2012._

This notebook is based on part of the tutorial "Optimised Symbolic Finite Difference Computation with Devito" presented at the Intel® HPC Developer Conference 2017.