---
title: "4.3.1 Implementation shortfall, VWAP, TWAP, Almgren-Chriss"
description: "Measure execution with implementation shortfall, compare VWAP and TWAP, estimate square-root impact and solve Almgren-Chriss with cvxpy."
url: "https://pythonparatrading.com/en/python-trading-course/level-4/implementation-shortfall-vwap-twap-almgren-chriss/"
lang: "en"
image: "https://pythonparatrading.com/og/en/python-trading-course/level-4/implementation-shortfall-vwap-twap-almgren-chriss.png"
---

# 4.3.1 · Implementation shortfall, VWAP, TWAP and Almgren-Chriss

> What it costs to execute a large order is not visible in the price of a single trade: it is spread across the bid-ask spread, the impact paid on the spot, the impact that stays in the price and the market move while the order executes. This lesson measures it on 727 orders of a fictitious desk with known impact, recovers the square-root law, checks why permanent impact is so hard to estimate and solves the Almgren-Chriss problem with the impact the data actually has.

- Course: [Level 4: Portfolios, risk, execution and live trading](https://pythonparatrading.com/en/python-trading-course/level-4/)
- Module: 4.3 · Microstructure and execution
- Duration: 70 min
- Prerequisites: [1.1.5 · Costs, settlement, leverage and margin](https://pythonparatrading.com/en/python-trading-course/level-1/trading-costs-settlement-leverage-and-margin/), [2.2.3 · Frequencies: intraday, swing and position trading, and the capacity of a strategy](https://pythonparatrading.com/en/python-trading-course/level-2/trading-frequency-and-strategy-capacity/), [2.4.7 · Cost modeling: commissions, bid-ask spread, impact, financing and securities borrowing](https://pythonparatrading.com/en/python-trading-course/level-2/backtest-transaction-cost-modeling/), [2.4.8 · Fill modeling: limit orders, queue, partial fills, gaps and stops](https://pythonparatrading.com/en/python-trading-course/level-2/backtest-order-fill-modeling/)

## By the end of this lesson

- Calculate Perold's implementation shortfall and break it down into delay, execution and opportunity cost.
- Compare an execution against arrival, VWAP, TWAP and the close, and explain why VWAP hides cost when the order is a large share of volume.
- Recover the square-root law (exponent and coefficient) from the child executions and see why the estimate with parent orders comes out biased.
- Estimate permanent impact and its error, and calculate how many orders would be needed to measure it precisely.
- Set up Almgren-Chriss, solve it with cvxpy with square-root temporary impact and compare it with the closed-form linear-impact formula.

## The paper portfolio and the real one

In 1988, André Perold proposed measuring execution by comparing two portfolios: the **paper** one, which buys or sells
everything at the price at the moment of the decision and with no costs, and the **real** one, which pays what it
pays and may not complete the order. The difference is the *implementation shortfall* (IS). Its parts tell a story
that the average execution price does not:

```mermaid
flowchart TB
  D[Decision] -- delay --> L[Arrival]
  L -- execution --> E[Average price]
  L -- opportunity --> C[Close: not done]
```

For an order of $Q$ shares with side $s$ (+1 buy, −1 sell), decision price $P_d$, market arrival price $P_0$,
execution of $q$ shares at an average price $\overline{P}$ and close $P_c$, as a fraction of the decision value:

$$
\text{IS} = \frac{s\,(P_0 - P_d)}{P_d}
+ \frac{s\,(\overline{P} - P_0)}{P_d}\,\frac{q}{Q}
+ \frac{s\,(P_c - P_0)}{P_d}\,\frac{Q - q}{Q}
$$

The three terms are, in order, delay, execution and opportunity. Delay is how much the price moves between the
decision and the sending of the order (the responsibility of whoever decides or of the process); execution, what it
costs to execute; opportunity, what is missed on what was not executed. Lesson 1.1.5 introduced the idea; here it is
measured.

The data is synthetic: four fictitious stocks (from `ARVA`, very liquid, to `DUNE`, illiquid and volatile) over 60
NYSE sessions of 2025, with 727 parent orders from a desk (TWAP, VWAP, POV and IS), their 92,409 child orders and the
market's 1-minute bars, which already carry the impact of those orders. Impact has known parameters
(`ejecucion4_acciones.csv`): temporary $\eta\,\sigma_d\,\rho^{0.5}$ with $\eta = 0.40$ and $\rho$ the order's
participation in the minute's volume, and linear permanent $\gamma\,\sigma_d\,q/\text{ADV}$ with $\gamma = 1.0$.

```python
import numpy as np
import pandas as pd

acc = pd.read_csv("ejecucion4_acciones.csv", index_col="ticker")
ord_ = pd.read_csv("ejecucion4_ordenes.csv",
                   parse_dates=["inicio", "fin_real"])
ord_["lado_s"] = np.where(ord_["lado"] == "compra", 1, -1)
ord_["sigma"] = ord_["ticker"].map(acc["sigma_diaria"])
ord_["adv"] = ord_["ticker"].map(acc["adv_acciones"])
ord_["minutos"] = (ord_["fin_real"] - ord_["inicio"]
                   ).dt.total_seconds() / 60
print(acc[["adv_acciones", "sigma_diaria"]])
ord_.groupby("algoritmo").agg(
    ordenes=("id_orden", "size"),
    pct_adv_mediana=("pct_adv", "median"),
    minutos_mediana=("minutos", "median"),
    completas=("cantidad_ejecutada",
               lambda q: (q == ord_.loc[q.index, "cantidad"]).mean()),
).round(3)
```
```output
        adv_acciones  sigma_diaria
ticker
ARVA        11000000         0.014
BELO         6500000         0.018
CIRA         1750000         0.023
DUNE          640000         0.030
           ordenes  pct_adv_mediana  minutos_mediana  completas
algoritmo
IS             154            0.010             75.5      1.000
POV            198            0.012             39.5      0.808
TWAP           171            0.012            106.0      1.000
VWAP           204            0.016             92.5      1.000
```

Orders range from 0.1% to 20% of average daily volume (ADV), with medians from 1% to 1.6% depending on the algorithm,
and last from half an hour to a whole session (medians from 40 to 106 minutes). TWAP, VWAP and IS orders always
complete; POV orders take a fixed share of volume and 19% do not finish before the close. The previous session's
close is used as the decision price (orders on the first day do not have one and are left out):

```python
def pb(x, ref):
    return 1e4 * x / ref


o = ord_.dropna(subset=["cierre_anterior"]).copy()
s, Q, q = o["lado_s"], o["cantidad"], o["cantidad_ejecutada"]
dec, lleg = o["cierre_anterior"], o["precio_llegada"]
o["retraso"] = pb(s * (lleg - dec), dec)
o["ejecucion"] = pb(s * (o["precio_medio"] - lleg), dec) * q / Q
o["oportunidad"] = pb(s * (o["cierre"] - lleg), dec) * (Q - q) / Q
o["is_total"] = o[["retraso", "ejecucion", "oportunidad"]].sum(axis=1)
print(len(o), "orders with a previous close")
o.groupby("algoritmo")[["retraso", "ejecucion", "oportunidad",
                        "is_total"]].mean().round(1)
```
```output
710 orders with a previous close
           retraso  ejecucion  oportunidad  is_total
algoritmo
IS             1.8       40.6          0.0      42.5
POV           17.2       33.1          2.4      52.7
TWAP          -8.0       47.2          0.0      39.2
VWAP          -2.5       44.8          0.0      42.3
```

The average IS is 39 to 53 bp depending on the algorithm, and almost all of it is execution (33 to 47 bp). Delay is
market noise between the close and the sending of the order: +17.2 bp in POV, −8.0 in TWAP, unrelated to the
algorithm. Opportunity only appears in POV (2.4 bp), the only one that leaves orders incomplete. It is a useful
reminder: an algorithm that completes less looks cheaper if you only look at what was executed.

## Four benchmark prices

The execution component can be measured against several benchmarks, and each one says something different:

| Benchmark | What it is | What it rewards |
|---|---|---|
| Arrival | Mid when the order is sent | Executing cheaply relative to the moment of the decision |
| Interval VWAP | Market volume-weighted average price while the order executes | Following the market |
| Interval TWAP | Simple average of per-minute prices over the interval | Spreading evenly over time |
| Close | Closing price | Matching whoever values the portfolio at the close |

```python
barras = pd.read_parquet("ejecucion4_barras_1min.parquet")
por_ticker = {t: g.set_index("ts") for t, g in barras.groupby("ticker")}


def twap_intervalo(r):
    b = por_ticker[r.ticker].loc[r.inicio:r.fin_real
                                 - pd.Timedelta(seconds=1)]
    return b["vwap"].mean()


ord_["twap_intervalo"] = ord_.apply(twap_intervalo, axis=1)
ref = {"llegada": "precio_llegada", "vwap": "vwap_intervalo",
       "twap": "twap_intervalo", "cierre": "cierre"}
for nombre, col in ref.items():
    ord_[f"frente_{nombre}"] = pb(ord_["lado_s"] * (
        ord_["precio_medio"] - ord_[col]), ord_[col])
cols = [f"frente_{n}" for n in ref]
print(ord_.groupby("algoritmo")[cols].mean().round(1))
ord_[cols].std().round(1)
```
```output
           frente_llegada  frente_vwap  frente_twap  frente_cierre
algoritmo
IS                   40.2         29.1         29.5            9.5
POV                  37.5         29.8         30.1           39.4
TWAP                 46.6         26.9         26.0           22.5
VWAP                 44.2         27.9         26.7           24.6
frente_llegada     72.6
frente_vwap        17.0
frente_twap        18.1
frente_cierre     140.3
dtype: float64
```

Against arrival, the orders cost 37.5 to 46.6 bp; against the interval VWAP or TWAP, 26.0 to 30.1: about 10 to 20 bp
less. The standard deviations explain why managers prefer VWAP: 17.0 bp against 72.6 for arrival, because VWAP moves
with the market during the order and removes almost all of its noise. The close is the noisiest benchmark
(140.3 bp): it says nothing about the quality of an execution that finished hours earlier.

What VWAP removes is not only noise:

```python
def volumen_intervalo(r):
    b = por_ticker[r.ticker].loc[r.inicio:r.fin_real
                                 - pd.Timedelta(seconds=1)]
    return b["volumen"].sum()


ord_["participacion_media"] = (ord_["cantidad_ejecutada"]
                               / ord_.apply(volumen_intervalo, axis=1))
tramos = pd.cut(ord_["participacion_media"],
                [0, 0.02, 0.05, 0.10, 0.20, 1.0])
ord_.groupby(tramos, observed=True)[
    ["frente_llegada", "frente_vwap"]].agg(["mean", "size"]).round(1)
```
```output
                    frente_llegada      frente_vwap
                              mean size        mean size
participacion_media
(0.0, 0.02]                   25.0  160        20.8  160
(0.02, 0.05]                  38.6  126        23.2  126
(0.05, 0.1]                   29.2  113        28.1  113
(0.1, 0.2]                    38.6  203        32.0  203
(0.2, 1.0]                    84.9  125        38.1  125
```

With less than 2% of volume, arrival and VWAP almost coincide (25.0 and 20.8 bp). With more than 20%, the order costs
84.9 bp against arrival and 38.1 against VWAP: half of the cost has disappeared from the measure. The order's own
executions are a large part of the interval's volume and their impact raises the VWAP; the benchmark moves with the
order. A trader evaluated against VWAP has no incentive to execute with less impact; one evaluated against arrival
does.

*Fig. — The desk's 727 orders by average share of volume: cost against the arrival price and against the interval VWAP. The more the order weighs in the volume, the more cost the VWAP hides.* Pairs of bars by participation band. Up to 2%: 25.0 bp against arrival and 20.8 against VWAP. 2-5%: 38.6 and 23.2. 5-10%: 29.2 and 28.1. 10-20%: 38.6 and 32.0. Over 20% (highlighted): 84.9 against arrival and 38.1 against VWAP.

## What is inside the cost

For each order, the generator stores the true decomposition of its shortfall against arrival
(`ejecucion4_ordenes_verdad.csv`): half spread, temporary impact, the permanent impact paid by its own child orders
and the mid move without impact (the market). In reality that decomposition is not known; here it shows what the
number is made of:

```python
verdad = pd.read_csv("ejecucion4_ordenes_verdad.csv")
ord_ = ord_.merge(verdad, on="id_orden")
partes = ["coste_horquilla_pb", "coste_temporal_pb",
          "coste_permanente_pb", "mercado_pb"]
tabla = ord_.groupby("ticker")[partes + ["frente_llegada"]].mean()
tabla.columns = ["horquilla", "temporal", "permanente", "mercado",
                 "IS"]
print(tabla.round(1))
ord_[partes + ["frente_llegada"]].std().round(1)
```
```output
        horquilla  temporal  permanente  mercado    IS
ticker
ARVA          0.3      18.0         1.9     -3.8  16.3
BELO          0.9      24.1         2.5      3.2  30.7
CIRA          4.8      32.0         3.5      2.8  43.0
DUNE         13.6      47.0         5.3     10.3  76.1
coste_horquilla_pb      5.5
coste_temporal_pb      18.7
coste_permanente_pb     4.7
mercado_pb             66.5
frente_llegada         72.6
dtype: float64
```

Cost grows from `ARVA` to `DUNE` in all three components: the spread (0.3 to 13.6 bp), the temporary impact (18.0 to
47.0) and the permanent impact (1.9 to 5.3). Temporary impact is the bulk of the cost in all four stocks. But the
standard deviations change the reading: 18.7 bp for the temporary impact, 4.7 for the permanent and 66.5 for the
market. Order by order, what decides the result is the market move during execution, which the trader does not
control. Measuring impact takes many orders.

Lesson 2.4.7 used the **square-root law** as an assumption, with its coefficient stated: the impact of an order of $Q$
shares executed over a day is approximately $Y\,\sigma_d \sqrt{Q/\text{ADV}}$, with $Y$ of the order of 1 in empirical
studies (Tóth et al., 2011, explain it and survey it). With the executions it can be estimated. Each child order pays,
against the mid prevailing just before, the half spread plus the minute's temporary impact,
$\eta\,\sigma_d\,\rho^{\beta}$; subtracting the half spread, a regression without a constant on $\sigma_d\sqrt{\rho}$
gives $\eta$, and the average costs by participation bucket on a log-log scale give the exponent $\beta$:

```python
hijas = pd.read_parquet("ejecucion4_hijas.parquet")
hijas = hijas.merge(ord_[["id_orden", "lado_s", "sigma"]],
                    on="id_orden")
hijas["minuto"] = hijas["ts"].dt.floor("min")
propias = hijas.groupby(["ticker", "minuto"])["cantidad"].sum()
vol = barras.set_index(["ticker", "ts"])["volumen"]
vol.index = vol.index.set_names(["ticker", "minuto"])
hijas = hijas.join(propias.rename("propias"), on=["ticker", "minuto"])
hijas = hijas.join(vol, on=["ticker", "minuto"])
hijas["rho"] = hijas["propias"] / hijas["volumen"]
hijas["coste"] = (hijas["lado_s"] * (hijas["precio"]
                                     - hijas["medio_antes"])
                  - hijas["horquilla"] / 2) / hijas["medio_antes"]

x = hijas["sigma"] * np.sqrt(hijas["rho"])
eta = (x * hijas["coste"]).sum() / (x**2).sum()
print(f"{len(hijas)} child orders; median participation "
      f"{hijas['rho'].median():.3f}; eta = {eta:.3f}")

hijas["tramo"] = pd.qcut(np.log(hijas["rho"]), 20)
cubos = hijas.groupby("tramo", observed=True).agg(
    rho=("rho", "mean"), coste_sigma=("coste", "sum"),
    n=("coste", "size"), sigma=("sigma", "sum"))
cubos["coste_sigma"] = cubos["coste_sigma"] / cubos["sigma"]
pendiente, corte = np.polyfit(np.log(cubos["rho"]),
                              np.log(cubos["coste_sigma"]), 1)
print(f"Exponent {pendiente:.3f}; coefficient {np.exp(corte):.3f}")
```
```output
92409 child orders; median participation 0.098; eta = 0.400
Exponent 0.500; coefficient 0.400
```

With 92,409 child orders the recovery is exact in practice: $\eta = 0.400$, exponent 0.500 and a log-fit coefficient
of 0.400, the generator's parameters. The figure at the top shows the twenty buckets on the line with slope 1/2. The
linear alternative (slope 1), the one the closed-form Almgren-Chriss formula assumes, only matches the square root at
one point: it exaggerates the cost of high participation and underestimates it at low participation.

Live, this estimate requires the data of each child order with the mid prevailing just before sending it, which a
broker does not always provide. With parent-order data:

```python
hs = hijas.assign(h=1e4 * hijas["horquilla"] / 2
                  / hijas["medio_antes"] * hijas["cantidad"])
hs = hs.groupby("id_orden")[["h", "cantidad"]].sum()
ord_ = ord_.join((hs["h"] / hs["cantidad"]).rename("horq_obs"),
                 on="id_orden")
ord_["impacto"] = ord_["frente_llegada"] - ord_["horq_obs"]


def sin_constante(y, x):
    b = (x * y).sum() / (x**2).sum()
    e = np.sqrt(((y - b * x)**2).sum() / (len(y) - 1)
                / (x**2).sum())
    return b, e


filas = {}
for alg, d in ord_.groupby("algoritmo"):
    x = 1e4 * d["sigma"] * np.sqrt(d["participacion_media"])
    b, e = sin_constante(d["impacto"], x)
    b2, _ = sin_constante(d["impacto"] - d["mercado_pb"], x)
    filas[alg] = {"Y": b, "error": e, "Y_sin_mercado": b2}
pd.DataFrame(filas).T.round(3)
```
```output
          Y  error  Y_sin_mercado
IS    0.570  0.046          0.538
POV   0.442  0.055          0.446
TWAP  0.584  0.087          0.497
VWAP  0.573  0.052          0.494
```

With parent orders, against arrival and subtracting the observable half spread, the coefficient comes out between
0.44 and 0.58 depending on the algorithm, with errors of 0.05 to 0.09. It is biased upward for two reasons: the parent
order also pays half of its own permanent impact, which grows with size, and market noise. Subtracting the true market
move (something that cannot be done in reality) the bias shrinks (0.45 to 0.54) but does not disappear, because the
permanent impact is still inside. POV comes out closer to 0.4 because its orders are short (40 minutes median) and
accumulate less permanent impact.

**Permanent** impact is what stays in the price after the order. With what is observable, it is estimated by
subtracting the spread and the temporary impact (the latter with the $\eta$ from the child orders) from the cost of
each order and comparing the remainder with $\sigma_d\,X/2$, half of the total permanent impact of an order of
$X = q/\text{ADV}$:

```python
t = hijas.assign(t=1e4 * eta * hijas["sigma"] * np.sqrt(hijas["rho"])
                 * hijas["cantidad"])
t = t.groupby("id_orden")["t"].sum() / hs["cantidad"]
ord_["temporal_est"] = ord_["id_orden"].map(t)
X = ord_["cantidad_ejecutada"] / ord_["adv"]
z = 1e4 * ord_["sigma"] * X / 2
resto = ord_["impacto"] - ord_["temporal_est"]
g, e = sin_constante(resto, z)
g2, e2 = sin_constante(resto - ord_["mercado_pb"], z)
print(f"gamma with what is observable: {g:.2f} ± {e:.2f}")
print(f"gamma subtracting the market (truth): {g2:.3f} ± {e2:.3f}")
n_necesarias = len(ord_) * (e / 0.10)**2
print(f"Orders for an error of ±0.10: {n_necesarias:,.0f}")
```
```output
gamma with what is observable: 1.76 ± 0.42
gamma subtracting the market (truth): 0.983 ± 0.008
Orders for an error of ±0.10: 12,881
```

With the 727 orders, $\gamma = 1.76 \pm 0.42$, when the truth is 1.0: within two standard errors, but useless for
deciding anything. Subtracting the true market move gives $0.983 \pm 0.008$: the permanent impact is there, covered by
noise tens of times larger. For an error of ±0.10 about 12,900 orders like these would be needed. That is why
estimates of permanent impact are made with hundreds of thousands of orders from many desks, and why in your own cost
model permanent impact is an assumption that has to be stated and varied, not a measured figure.

<aside class="callout" data-kind="aviso">

An impact coefficient estimated with few orders from a single desk is not a property of the market: it depends on the
algorithms, sizes and dates of that sample. Before using it to size a strategy (lesson 2.2.3), compare it with the
range in the literature and test the result with double and half its value.

</aside>

## Urgency: the Almgren and Chriss problem

Executing fast raises temporary impact (more participation); executing slowly leaves the order exposed to the market
move for longer. Almgren and Chriss (2000) framed it as a portfolio problem: choose the execution schedule $x_k$ (what
remains to be executed at the end of each slice) that minimizes

$$
E[\text{cost}] + \lambda\,\text{Var}[\text{cost}], \qquad
\text{Var}[\text{cost}] = \sigma^2 \sum_k \tau\, x_k^2
$$

with $\lambda$ the risk aversion. With temporary impact **linear** in the execution speed and constant liquidity, the
solution is closed-form: $x(t) = X \sinh(\kappa (T - t)) / \sinh(\kappa T)$, with $\kappa^2 = \lambda \sigma^2 / \tilde\eta$.
With $\lambda = 0$ it gives TWAP; the larger $\lambda$, the more execution is front-loaded. The desk's IS orders follow
that trajectory with urgencies $\kappa \in \{0.5,\ 2,\ 5\}$. With the true decomposition the trade-off is visible:

```python
d = ord_[ord_["algoritmo"] == "IS"].copy()
d["mercado_en_sigmas"] = d["mercado_pb"] / (1e4 * d["sigma"])
d.groupby("urgencia_kappa").agg(
    ordenes=("id_orden", "size"), minutos=("minutos", "median"),
    temporal=("coste_temporal_pb", "mean"),
    riesgo=("mercado_en_sigmas", "std")).round(3)
```
```output
                ordenes  minutos  temporal  riesgo
urgencia_kappa
0.5                  40     95.5    25.028   0.237
2.0                  53     98.0    31.611   0.233
5.0                  61     60.0    38.690   0.136
```

Going from $\kappa = 0.5$ to 5, the average temporary impact rises from 25.0 to 38.7 bp and the risk, measured as the
standard deviation of the market move during the order in daily standard deviations, falls from 0.237 to 0.136.
Urgent orders last less (60 minutes median against 96) and pay more for participating more. Between 0.5 and 2 the
risk barely changes: the two lowest urgencies have similar durations.

The data does not meet the assumption of the closed-form formula: temporary impact is $\eta\,\sigma_d\,\rho^{0.5}$, in
the square root of participation, and volume follows a U-shaped curve. The cost of a slice $k$ with $n_k$ shares and
a market volume $V_k$ is $\eta\,\sigma_d\, n_k^{1.5} / V_k^{0.5}$: convex in $n_k$, so the problem is still convex and
cvxpy solves it with `cp.power(u, 1.5)` (a power cone that CLARABEL supports). The case: sell 10% of `DUNE`'s ADV in
one session, in 13 slices of 30 minutes, with the average volume curve of the data, the generator's $\eta$ and
$\gamma$ and the Almgren and Chriss variance with volatility spread evenly over the session (the data has it U-shaped:
it is a simplification). With linear permanent impact, the permanent cost does not depend on the schedule: it always
adds up to $\gamma\,\sigma_d\,X/2$.

```python
import cvxpy as cp

T_ = "DUNE"
b = por_ticker[T_]
minuto = (b.index.hour * 60 + b.index.minute) - (13 * 60 + 30)
curva = b.groupby([b.index.date, minuto // 30])["volumen"].sum()
V = curva.groupby(level=1).mean().to_numpy()   # 13 slices of 30 min
N, tau = len(V), 1 / len(V)                    # tau in days
ADV, sig = acc.at[T_, "adv_acciones"], acc.at[T_, "sigma_diaria"]
Xq = 0.10 * ADV                                # 10% of ADV
ETA, GAMMA = 0.40, 1.0                         # the generator's


def coste(u):
    """Expected cost and deviation (bp) of schedule u (fraction of
    the order in each slice). Temporary impact in the square root of
    participation, linear permanent, variance of what remains."""
    temporal = 1e4 * ETA * sig * np.sum(u**1.5 * np.sqrt(Xq / V))
    permanente = 1e4 * GAMMA * sig * Xq / ADV / 2
    resto = 1 - np.cumsum(u)
    desv = 1e4 * sig * np.sqrt(tau * np.sum(resto**2))
    return temporal + permanente, desv


def optimo(lam):
    u = cp.Variable(N, nonneg=True)
    temporal = 1e4 * ETA * sig * cp.sum(
        cp.multiply(np.sqrt(Xq / V), cp.power(u, 1.5)))
    varianza = (1e4 * sig)**2 * tau * cp.sum_squares(1 - cp.cumsum(u))
    cp.Problem(cp.Minimize(temporal + lam * varianza),
               [cp.sum(u) == 1]).solve(solver="CLARABEL")
    return np.clip(u.value, 0, None)


print("Average volume per slice (thousands):", np.round(V / 1e3, 1))
cal = {"TWAP": np.full(N, 1 / N), "VWAP": V / V.sum(),
       "optimo_lam0": optimo(0.0)}
pd.DataFrame({k: coste(u) for k, u in cal.items()},
             index=["coste_pb", "desv_pb"]).T.round(2)
```
```output
Average volume per slice (thousands): [98.9 68.2 51.2 45.5 44.3 43.1 44.  43.4 44.9 44.5 47.4 52.5 92. ]
             coste_pb  desv_pb
TWAP            51.89   163.18
VWAP            50.78   153.48
optimo_lam0     50.78   153.48
```

VWAP costs 50.8 bp and TWAP 51.9: with square-root impact, the optimum without risk aversion is **VWAP**, not TWAP.
Total cost is proportional to $\sum_k n_k^{1.5} / V_k^{0.5}$, and minimizing it with $\sum n_k$ fixed gives
$n_k \propto V_k$: constant participation. The numerical solution with $\lambda = 0$ confirms it (the same 50.78 bp).
The standard deviation of the cost, 153 bp, is three times the expected cost: an order of 10% of the ADV of a stock
with 3% daily volatility is, above all, an involuntary bet on the price of the session.

To use the closed-form formula, impact has to be **linearized**: choose a linear $\tilde\eta$ that gives the same cost
per share as the square root with the TWAP schedule, and apply the hyperbolic sine formula. The two schedules are
compared for several values of $\lambda$, both evaluated with the data's square-root impact:

```python
# linear eta that gives the same cost per share as the square root with TWAP
n_tw = Xq / N
eta_lin = ETA * sig * tau / np.sqrt(n_tw * V.mean())


def cerrado(lam):
    """Almgren-Chriss (2000) schedule with linear impact:
    x(t) = sinh(k (T - t)) / sinh(k T), T = 1 day."""
    k = np.sqrt(lam * 1e4 * sig**2 / (eta_lin * Xq))
    t = np.arange(N + 1) / N
    x = np.sinh(k * (1 - t)) / np.sinh(k) if k > 1e-9 else 1 - t
    return -np.diff(x)


filas = {}
for lam in (1e-4, 1e-3, 3e-3, 1e-2, 3e-2):
    for nombre, u in (("raiz_cvxpy", optimo(lam)),
                      ("lineal_cerrado", cerrado(lam))):
        c, sd = coste(u)
        filas[(lam, nombre)] = {"coste_pb": c, "desv_pb": sd,
                                "objetivo": c + lam * sd**2,
                                "primer_tramo": u[0]}
frontera = pd.DataFrame(filas).T
frontera.round(3)
```
```output
                       coste_pb  desv_pb  objetivo  primer_tramo
0.0001 raiz_cvxpy        50.911  144.054    52.986         0.160
       lineal_cerrado    51.809  160.198    54.375         0.083
0.0010 raiz_cvxpy        55.999   96.594    65.329         0.306
       lineal_cerrado    52.101  139.021    71.428         0.126
0.0030 raiz_cvxpy        65.556   61.825    77.024         0.465
       lineal_cerrado    54.772  112.042    92.432         0.192
0.0100 raiz_cvxpy        80.118   32.253    90.520         0.669
       lineal_cerrado    62.446   77.098   121.888         0.320
0.0300 raiz_cvxpy        93.072   15.553   100.329         0.826
       lineal_cerrado    72.325   49.674   146.349         0.487
```

With the same $\lambda$, the schedule from the linear formula is much less urgent than the square-root optimum: with
$\lambda = 0.01$ it executes 32% in the first half hour against 67%, and its objective (cost plus $\lambda$ times the
variance) is 121.9 against 90.5. The reason lies in the shape of the impact: with linear impact, the total cost of a
slice grows with $n^2$; with the square root, with $n^{1.5}$. The linear formula exaggerates what it costs to
concentrate execution and delays it.

*Fig. — Schedules to sell 10% of DUNE's ADV in one session: VWAP, TWAP and the Almgren-Chriss optimum with aversion λ = 0.01, solved with square-root impact (cvxpy) and with the closed-form linear-impact formula.* Four curves of the fraction of the order still to execute, from 9:30 to 16:00 in 30-minute slots. TWAP (dashed) falls in a straight line. VWAP falls faster in the first and last half hour, where volume is higher. The square-root optimum with λ = 0.01 (highlighted) executes 67% in the first half hour and almost everything before noon. The linear closed form with the same λ (blue) is much slower: 32% in the first half hour.

As a family, on the other hand, the hyperbolic sine schedules are not far from the frontier: they sit between 1 and
2 bp above it for the same risk. The linearization error is not so much the shape of the schedule as the urgency it
assigns to each $\lambda$: whoever calibrates their risk aversion with the linear formula and executes with
square-root impact gets a slower and riskier schedule than intended.

*Fig. — Expected cost against risk of the schedules to sell 10% of DUNE's ADV, evaluated with the data's square-root impact: the numerical frontier and the closed-form linear schedules.* Two curves on the standard deviation (horizontal) and expected cost (vertical) plane, for 18 risk aversions. The square-root frontier (highlighted) runs from VWAP, 50.8 bp with 153 of deviation, to about 100 bp with less than 10 of deviation. The linear-formula schedules (blue) lie above and to the right: for the same risk they cost more. TWAP lies to the right of VWAP: 51.9 bp with 163.

## Common mistakes

- **Measuring only what was executed.** Without the opportunity cost, an algorithm that does not complete looks cheap.
- **Evaluating large orders against the interval VWAP**: the benchmark includes their own impact.
- **Confusing noise with cost.** With 73 bp of standard deviation per order, comparing two algorithms with a few dozen
  orders says nothing; you need the average of many and its standard error.
- **Estimating impact with parent orders** without separating the permanent impact: the coefficient comes out
  inflated.
- **Treating permanent impact as measured.** With a single desk it is an assumption.
- **Using the closed-form Almgren-Chriss formula** with square-root impact without saying so: the urgency comes out
  wrong.
- **Forgetting the volume curve.** With U-shaped liquidity, the risk-neutral schedule is VWAP, not TWAP.

## Exercises

1. Estimate the coefficient $Y$ of the square-root law in its classic form, $Y\,\sigma_d\sqrt{Q/\text{ADV}}$, using
   only the VWAP orders that last the whole session (390 minutes or almost). How many are there and what is the error
   of the estimate?

**Solution.**


```python
dia = ord_[(ord_["algoritmo"] == "VWAP") & (ord_["minutos"] >= 380)]
x = 1e4 * dia["sigma"] * np.sqrt(dia["cantidad_ejecutada"]
                                 / dia["adv"])
Y, err = sin_constante(dia["impacto"], x)
Y2, _ = sin_constante(dia["impacto"] - dia["mercado_pb"], x)
print(f"{len(dia)} one-day orders: Y = {Y:.2f} ± {err:.2f};"
      f" without the market, {Y2:.2f}")
```
```output
18 one-day orders: Y = 0.94 ± 0.81; without the market, 0.59
```

There are only 18 VWAP orders that last the whole session, and the estimate is $Y = 0.94 \pm 0.81$: the error is
almost as large as the coefficient. Subtracting the true market move gives 0.59, above 0.40 because the parent order
also pays half of its permanent impact. The classic form of the law, the one used to size a strategy, cannot be
estimated with one desk: that is why the child orders, with their minute-by-minute participation, are the good
source.

2. Repeat the Almgren-Chriss optimum with $\lambda = 0.01$ for 10% of the ADV of `ARVA`, the most liquid and least
   volatile stock. How much does it execute in the first half hour and how much does it cost compared with `DUNE`?

**Solution.**


```python
def caso(ticker, pct=0.10, lam=1e-2):
    global V, ADV, sig, Xq
    b = por_ticker[ticker]
    m = (b.index.hour * 60 + b.index.minute) - (13 * 60 + 30)
    curva = b.groupby([b.index.date, m // 30])["volumen"].sum()
    V = curva.groupby(level=1).mean().to_numpy()
    ADV = acc.at[ticker, "adv_acciones"]
    sig = acc.at[ticker, "sigma_diaria"]
    Xq = pct * ADV
    u = optimo(lam)
    c, sd = coste(u)
    return {"primer_tramo": u[0], "coste_pb": c, "desv_pb": sd,
            "vwap_pb": coste(V / V.sum())[0]}


pd.DataFrame({t: caso(t) for t in ("ARVA", "DUNE")}).T.round(3)
```
```output
      primer_tramo  coste_pb  desv_pb  vwap_pb
ARVA         0.547    33.311   23.138   23.912
DUNE         0.669    80.118   32.253   50.776
```

With the same risk aversion, `ARVA` executes 54.7% in the first half hour against 66.9% for `DUNE`, and it costs
33.3 bp with a deviation of 23.1, against 80.1 and 32.3. Its VWAP would cost 23.9 bp (50.8 for `DUNE`). With less
volatility, the risk of waiting weighs less and the schedule is front-loaded less; with $\lambda$ in bp⁻¹, the same
aversion means different things for different stocks, and in practice it is set per order or by the money size of the
risk that is accepted.

3. A manager decides to buy at yesterday's close and the trader receives the order at noon. The desk's average delay
   is several basis points, positive or negative depending on the day. Who should be accountable for that cost and how
   would you measure it to separate it from execution?

**Solution.**

The delay between the decision and arrival does not depend on the trader, but on the process: when the decision is
made, when it is approved and when it is sent. It is measured separately with Perold's IS, which separates it from
execution, and is evaluated on average over many orders: order by order it is market noise (on the desk, +17.2 bp in
POV and −8.0 in TWAP, unrelated to the algorithm). If on average it is positive and significant on buys and negative
on sells, the price moves against you between the decision and the sending of the order: information is leaking or a
signal that others see earlier, and the solution is to shorten the process, not to change the algorithm. The trader
is accountable for execution against arrival; the manager, for the delay and for the opportunity cost of what they
decided not to execute.

## Self-assessment

1. A buy order is decided at yesterday's close (100.00), reaches the market at 100.20 and half of it executes at an average of 100.50; the rest does not execute and the stock closes at 101.00. What is the implementation shortfall in basis points over the decision price?

   - A. 50 bp: only what was executed counts
   - B. 75 bp: 20 of delay, 15 of execution and 40 of opportunity
   - C. 30 bp: execution against arrival
   - D. 100 bp: from yesterday's close to today's

   *Answer B: 75 bp: 20 of delay, 15 of execution and 40 of opportunity.* Delay: 0.20 over 100 = 20 bp on the whole order. Execution: (100.50 − 100.20) on the executed half = 30 bp × 0.5 = 15. Opportunity: (101.00 − 100.20) on the unexecuted half = 80 × 0.5 = 40. Total, 75 bp: how much the real portfolio underperforms the paper portfolio.

2. In the data, orders that are more than 20% of the volume of their interval cost 84.9 bp against arrival and 38.1 against the interval VWAP. Why such a difference?

   - A. Because VWAP is a bad calculation
   - B. Because those orders execute at the close
   - C. Because the bid-ask spread is wider
   - D. Because the order's own executions are part of the VWAP: the order's impact raises the market's average price and the benchmark moves with it

   *Answer D: Because the order's own executions are part of the VWAP: the order's impact raises the market's average price and the benchmark moves with it.* With high participation, a good part of the interval's volume is the order itself. The interval VWAP already includes its impact, so comparing against it rewards whoever moves the price. Arrival does not depend on the order.

3. The regression on the 92,409 child orders gives a temporary impact of 0.400·σ·ρ^0.500. With parent orders, the same coefficient comes out between 0.44 and 0.58. Why?

   - A. Because the child orders have errors
   - B. Because parent orders do not follow the square-root law
   - C. Because the cost of a parent order also includes the permanent impact accumulated during the order and market noise, which the regression partly attributes to the temporary one
   - D. Because of rounding to the cent

   *Answer C: Because the cost of a parent order also includes the permanent impact accumulated during the order and market noise, which the regression partly attributes to the temporary one.* Each child order is measured against the mid prevailing just before, which already includes the accumulated permanent impact; the parent order is measured against arrival, and also pays half of its permanent impact and the market move. Subtracting the market move (which is not known in reality) the coefficient drops to 0.45-0.54.

4. With what is observable, the permanent impact of the 727 orders comes out at γ = 1.76 ± 0.42 (truth: 1.0). How many orders would be needed for an error of ±0.10?

   - A. About 13,000: the error falls with the square root of the number of orders
   - B. About 1,000
   - C. 727 are enough if the child orders are used
   - D. It is impossible to measure

   *Answer A: About 13,000: the error falls with the square root of the number of orders.* An error of 0.42 with 727 orders; for 0.10 it must be divided by 4.2, which takes 4.2² ≈ 17.7 times more orders: about 12,900. That is why impact studies use hundreds of thousands of orders from many desks.

5. With temporary impact proportional to the square root of participation and no risk aversion, which schedule minimizes expected cost?

   - A. VWAP: a constant participation in volume, because the cost per share only depends on participation
   - B. TWAP: the same number of shares every minute
   - C. Everything at the start
   - D. Everything at the close

   *Answer A: VWAP: a constant participation in volume, because the cost per share only depends on participation.* Total cost is proportional to the sum of n^1.5/V^0.5. Minimizing it with the sum of n fixed gives n proportional to V: constant participation. The Almgren-Chriss formula with linear impact and constant liquidity gives TWAP; with the U-shaped volume curve of the data, VWAP costs 50.8 bp against 51.9 for TWAP.

## Key points

- Perold's implementation shortfall compares the real portfolio with the paper one at decision prices and breaks down into delay, execution and opportunity. On the lesson's desk, about 42 bp on average with a standard deviation of 73: market noise dominates order by order.
- The interval VWAP is a convenient but soft benchmark: it includes the order's own executions. Above 20% of volume, orders cost 84.9 bp against arrival and only 38.1 against VWAP.
- The square-root law is recovered precisely from the child orders: exponent 0.500 and coefficient 0.400, those of the generator. With parent orders the coefficient comes out between 0.44 and 0.58, because it mixes with permanent impact and noise.
- Permanent impact is the hard parameter: 1.76 ± 0.42 with what is observable against 1.0 in truth; only by subtracting the market move, which is not known live, does it come out at 0.98. Measuring it with an error of ±0.10 would take about 13,000 orders.
- Almgren-Chriss balances expected cost and risk. With square-root temporary impact there is no closed form: it is solved with cvxpy (`cp.power(u, 1.5)`). Without risk aversion the optimum is VWAP; with aversion, it front-loads execution much more than the linear formula with the same λ.

## Resources

### Research

- [Quantitative Models and Algorithmic Execution](https://pythonparatrading.com/en/quantitative-models-algorithmic-execution/) · Research-0009 · 11 min

### Documentation

- [cvxpy: atomic functions (power, sum_squares, cumsum)](https://www.cvxpy.org/tutorial/functions/index.html) · cvxpy.org

### Books

- [A. F. Perold, “The implementation shortfall: Paper versus reality”, Journal of Portfolio Management 14(3) (1988), 4-9 (DOI)](https://doi.org/10.3905/jpm.1988.409150) · doi.org
- [R. Almgren and N. Chriss, “Optimal execution of portfolio transactions”, Journal of Risk 3(2) (2000/2001), 5-39 (DOI)](https://doi.org/10.21314/JOR.2001.041) · doi.org
- [B. Tóth et al., “Anomalous price impact and the critical nature of liquidity in financial markets”, Physical Review X 1 (2011) (DOI)](https://doi.org/10.1103/PhysRevX.1.021006) · doi.org

### Data

- [Level 4 data: the desk's parent orders (CSV)](https://pythonparatrading.com/curso/datos/nivel-4/ejecucion4_ordenes.csv) · Download · ejecucion4_ordenes.csv
- [Level 4 data: child executions (Parquet)](https://pythonparatrading.com/curso/datos/nivel-4/ejecucion4_hijas.parquet) · Download · ejecucion4_hijas.parquet

Previous lesson: [4.2.6 · Derivatives risk: portfolio Greeks, margin and liquidations](https://pythonparatrading.com/en/python-trading-course/level-4/derivatives-risk-greeks-margin-and-liquidations/)

Next lesson: [4.3.2 · Order flow imbalance, microstructure signals and latency](https://pythonparatrading.com/en/python-trading-course/level-4/order-flow-imbalance-and-latency/)