"""
## GAMSSOURCE: https://www.gams.com/latest/gamslib_ml/libhtml/gamslib_dyncge.html
## LICENSETYPE: Demo
## MODELTYPE: NLP
## KEYWORDS: nonlinear programming, general equilibrium model, social accounting matrix
A Recursive-Dynamic Standard CGE Model (DYNCGE)
This model is featured in the following book.
Hosoe, N., Gasawa, K., Hashimoto, H. Textbook of Computable General
Equilibrium Modeling: Programming and Simulations, 2nd Edition,
University of Tokyo Press. (in Japanese)
"""
from __future__ import annotations
import os
import numpy as np
from gamspy import (
Alias,
Container,
Equation,
Model,
Ord,
Parameter,
Problem,
Product,
Sense,
Set,
Sum,
Variable,
)
def main():
m = Container(
system_directory=os.getenv("SYSTEM_DIRECTORY", None),
)
# ===============================================================
# Definition of sets for suffix ---------------------------------
# ===============================================================
# Sets
u = Set(
m,
name="u",
records=[
"AGR",
"LMN",
"HMN",
"SRV",
"CAP",
"LAB",
"HOH",
"GOV",
"INV",
"EXT",
"IDT",
"TRF",
],
description="SAM entry",
)
i = Set(
m,
name="i",
domain=u,
records=["AGR", "LMN", "HMN", "SRV"],
description="goods",
)
h = Set(
m, name="h", domain=u, records=["CAP", "LAB"], description="factor"
)
h_mob = Set(
m, name="h_mob", domain=h, records=["LAB"], description="mobile factor"
)
t = Set(
m, name="t", records=[str(i) for i in range(31)], description="time"
)
# Alias
v = Alias(m, name="v", alias_with=u)
j = Alias(m, name="j", alias_with=i)
k = Alias(m, name="k", alias_with=h)
# ===============================================================
# Data for Dynamics ---------------------------------------------
# ===============================================================
# Scalar
ror = Parameter(m, name="ror", description="rate of return of capital")
dep = Parameter(m, name="dep", description="depreciation rate")
pop = Parameter(m, name="pop", description="population growth rate")
zeta = Parameter(
m,
name="zeta",
description="elasticity parameter for investment allocation",
)
ror[...] = 0.05
dep[...] = 0.04
pop[...] = 0.02
zeta[...] = 1
sam_data = np.array(
[
[
1643.017,
7560.896,
237.841,
1409.202,
0,
0,
3563.257,
0,
919.745,
62.464,
0,
0,
],
[
1485.854,
10803.527,
15330.764,
18597.270,
0,
0,
32220.169,
329.469,
802.026,
1196.525,
0,
0,
],
[
1071.954,
4277.721,
113390.269,
48734.424,
0,
0,
27648.678,
4.931,
34979.803,
55083.516,
0,
0,
],
[
2002.380,
11406.260,
50513.476,
177675.714,
0,
0,
234243.865,
90707.177,
79169.426,
17426.156,
0,
0,
],
[
5082.506,
7042.697,
21058.821,
163045.396,
0,
0,
0,
0,
0,
0,
0,
0,
],
[
1435.010,
8942.365,
42510.123,
222732.700,
0,
0,
0,
0,
0,
0,
0,
0,
],
[0, 0, 0, 0, 196229.420, 275620.198, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 52243.041, 0, 0, 0, 34024.445, 4774.091],
[0, 0, 0, 0, 0, 0, 121930.608, 0, 0, -6059.608, 0, 0],
[
2092.569,
23796.669,
30982.559,
10837.256,
0,
0,
0,
0,
0,
0,
0,
0,
],
[433.854, 4068.616, 9418.058, 20103.917, 0, 0, 0, 0, 0, 0, 0, 0],
[149.278, 2866.853, 1749.385, 8.575, 0, 0, 0, 0, 0, 0, 0, 0],
]
)
# ===============================================================
# SAM Data
# ===============================================================
SAM = Parameter(
m,
name="SAM",
domain=[u, v],
records=sam_data,
description="social accounting matrix for 2005 [bil. JPY]",
)
# Source: compiled by N. Hosoe, based on the I/O table for 2005
SAMGAP = Parameter(
m,
name="SAMGAP",
domain=u,
description="gaps between row sums and column sums",
)
SAMGAP[u] = Sum(v, SAM[u, v] - SAM[v, u])
# print(SAMGAP.records)
# ===============================================================
# Loading the initial values ------------------------------------
# ===============================================================
# Base year values
Y00 = Parameter(m, name="Y00", domain=j, description="composite factor")
F00 = Parameter(m, name="F00", domain=[h, j], description="factor input")
X00 = Parameter(
m, name="X00", domain=[i, j], description="intermediate input"
)
Z00 = Parameter(m, name="Z00", domain=j, description="gross output")
Xp00 = Parameter(
m, name="Xp00", domain=i, description="household consumption"
)
Xg00 = Parameter(
m, name="Xg00", domain=i, description="government consumption"
)
Xv00 = Parameter(m, name="Xv00", domain=i, description="investment demand")
E00 = Parameter(m, name="E00", domain=i, description="exports")
M00 = Parameter(m, name="M00", domain=i, description="imports")
Q00 = Parameter(
m, name="Q00", domain=i, description="Armington's composite good"
)
D00 = Parameter(m, name="D00", domain=i, description="domestic good")
Sp00 = Parameter(m, name="Sp00", description="private savings")
Td00 = Parameter(m, name="Td00", description="direct tax")
Tz00 = Parameter(m, name="Tz00", domain=j, description="production tax")
Tm00 = Parameter(m, name="Tm00", domain=j, description="import tariff")
III00 = Parameter(m, name="III00", description="composite investment")
II00 = Parameter(
m, name="II00", domain=j, description="sectoral investment"
)
KK00 = Parameter(m, name="KK00", domain=j, description="capital stock")
CC00 = Parameter(
m, name="CC00", description="composite consumption or felicity"
)
FF00 = Parameter(m, name="FF00", domain=h, description="factor endowment")
Sf00 = Parameter(
m, name="Sf00", description="foreign savings in US dollars"
)
tauz00 = Parameter(
m, name="tauz00", domain=i, description="production tax rate"
)
taum00 = Parameter(
m, name="taum00", domain=i, description="import tariff rate"
)
# Base run value
Y0 = Parameter(m, name="Y0", domain=[j, t], description="composite factor")
F0 = Parameter(m, name="F0", domain=[h, j, t], description="factor input")
X0 = Parameter(
m, name="X0", domain=[i, j, t], description="intermediate input"
)
Z0 = Parameter(m, name="Z0", domain=[j, t], description="gross output")
Xp0 = Parameter(
m, name="Xp0", domain=[i, t], description="household consumption"
)
Xv0 = Parameter(
m, name="Xv0", domain=[i, t], description="investment demand"
)
E0 = Parameter(m, name="E0", domain=[i, t], description="exports")
M0 = Parameter(m, name="M0", domain=[i, t], description="imports")
Q0 = Parameter(
m, name="Q0", domain=[i, t], description="Armington's composite good"
)
D0 = Parameter(m, name="D0", domain=[i, t], description="domestic good")
Sp0 = Parameter(m, name="Sp0", domain=t, description="private savings")
Td0 = Parameter(m, name="Td0", domain=t, description="direct tax")
Tz0 = Parameter(m, name="Tz0", domain=[j, t], description="production tax")
Tm0 = Parameter(m, name="Tm0", domain=[j, t], description="import tariff")
III0 = Parameter(
m, name="III0", domain=t, description="composite investment"
)
II0 = Parameter(
m, name="II0", domain=[j, t], description="sectoral investment"
)
KK0 = Parameter(m, name="KK0", domain=[j, t], description="capital stock")
CC0 = Parameter(
m,
name="CC0",
domain=t,
description="composite consumption or felicity",
)
FF0 = Parameter(
m, name="FF0", domain=[h, t], description="factor endowment"
)
pf0 = Parameter(
m, name="pf0", domain=[h, j, t], description="factor price"
)
py0 = Parameter(
m, name="py0", domain=[j, t], description="composite factor price"
)
pz0 = Parameter(
m, name="pz0", domain=[j, t], description="gross output price"
)
pq0 = Parameter(
m,
name="pq0",
domain=[i, t],
description="Armington's composite good price",
)
pe0 = Parameter(
m,
name="pe0",
domain=[i, t],
description="export price in local currency",
)
pm0 = Parameter(
m,
name="pm0",
domain=[i, t],
description="import price in local currency",
)
pd0 = Parameter(
m, name="pd0", domain=[i, t], description="domestic good price"
)
pk0 = Parameter(
m, name="pk0", domain=t, description="composite investment goods price"
)
epsilon0 = Parameter(
m, name="epsilon0", domain=t, description="exchange rate"
)
PRICE0 = Parameter(
m, name="PRICE0", domain=t, description="numeraire price"
)
# Exogenous variables
Xg0 = Parameter(
m, name="Xg0", domain=[i, t], description="government consumption"
)
Sf0 = Parameter(
m, name="Sf0", domain=t, description="foreign savings in US dollars"
)
pWe = Parameter(
m, name="pWe", domain=i, description="export price in US dollars"
)
pWm = Parameter(
m, name="pWm", domain=i, description="import price in US dollars"
)
tauz = Parameter(
m, name="tauz", domain=i, description="production tax rate"
)
taum = Parameter(
m, name="taum", domain=i, description="import tariff rate"
)
# for result reporting
Y1 = Parameter(m, name="Y1", domain=[j, t], description="composite factor")
F1 = Parameter(m, name="F1", domain=[h, j, t], description="factor input")
X1 = Parameter(
m, name="X1", domain=[i, j, t], description="intermediate input"
)
Z1 = Parameter(m, name="Z1", domain=[j, t], description="gross output")
Xp1 = Parameter(
m, name="Xp1", domain=[i, t], description="household consumption"
)
Xv1 = Parameter(
m, name="Xv1", domain=[i, t], description="investment demand"
)
E1 = Parameter(m, name="E1", domain=[i, t], description="exports")
M1 = Parameter(m, name="M1", domain=[i, t], description="imports")
Q1 = Parameter(
m, name="Q1", domain=[i, t], description="Armington's composite good"
)
D1 = Parameter(m, name="D1", domain=[i, t], description="domestic good")
Sp1 = Parameter(m, name="Sp1", domain=t, description="private saving")
Td1 = Parameter(m, name="Td1", domain=t, description="direct tax")
Tz1 = Parameter(m, name="Tz1", domain=[j, t], description="production tax")
Tm1 = Parameter(m, name="Tm1", domain=[i, t], description="import tariff")
FF1 = Parameter(
m,
name="FF1",
domain=[h, t],
description="initial sectoral factor uses",
)
II1 = Parameter(
m, name="II1", domain=[j, t], description="sectoral investment"
)
III1 = Parameter(
m, name="III1", domain=t, description="composite investment"
)
KK1 = Parameter(
m, name="KK1", domain=[j, t], description="sectoral capital stock"
)
CC1 = Parameter(m, name="CC1", domain=t, description="utility")
tauz1 = Parameter(
m, name="tauz1", domain=[i, t], description="production tax rates"
)
taum1 = Parameter(
m, name="taum1", domain=[i, t], description="import tariff rates"
)
pz1 = Parameter(
m, name="pz1", domain=[j, t], description="gross output price"
)
pd1 = Parameter(
m, name="pd1", domain=[j, t], description="domestic good price"
)
pm1 = Parameter(m, name="pm1", domain=[j, t], description="import price")
pe1 = Parameter(m, name="pe1", domain=[j, t], description="export price")
pq1 = Parameter(
m,
name="pq1",
domain=[j, t],
description="Armington's composite good price",
)
pf1 = Parameter(
m, name="pf1", domain=[h, j, t], description="factor price"
)
py1 = Parameter(
m, name="py1", domain=[j, t], description="composite factor price"
)
epsilon1 = Parameter(
m, name="epsilon1", domain=t, description="foreign exchange rate"
)
pk1 = Parameter(m, name="pk1", domain=t, description="capital good price")
PRICE1 = Parameter(
m, name="PRICE1", domain=t, description="numeraire price"
)
Td00[...] = SAM["GOV", "HOH"]
Tz00[j] = SAM["IDT", j]
Tm00[j] = SAM["TRF", j]
F00[h, j] = SAM[h, j]
Y00[j] = Sum(h, F00[h, j])
X00[i, j] = SAM[i, j]
Z00[j] = Y00[j] + Sum(i, X00[i, j])
M00[i] = SAM["EXT", i]
tauz00[j] = Tz00[j] / Z00[j]
taum00[j] = Tm00[j] / M00[j]
Xp00[i] = SAM[i, "HOH"]
CC00[...] = Sum(i, Xp00[i])
FF00[h] = SAM["HOH", h]
E00[i] = SAM[i, "EXT"]
D00[i] = (1 + tauz00[i]) * Z00[i] - E00[i]
Q00[i] = (1 + taum00[i]) * M00[i] + D00[i]
Sf00[...] = SAM["INV", "EXT"]
# ===============================================================
# Adjusting Investment in the SAM for the Assumed BAU Growth Path
# ===============================================================
# Scalars
III_ASS = Parameter(
m,
name="III_ASS",
description="required investment for the assumed growth",
)
III_SAM = Parameter(
m, name="III_SAM", description="observed investment in the SAM"
)
adj = Parameter(
m, name="adj", description="III_ASS vs. III_SAM [>1:more than actual]"
)
III_ASS[...] = (pop + dep) / ror * FF00["CAP"]
III_SAM[...] = Sum(i, SAM[i, "INV"])
adj[...] = III_ASS / III_SAM
# Adjusting investment level
Xv00[i] = SAM[i, "INV"] * adj
# Reallocating the gap made by the inv. adjustment to gov. cons.
Xg00[i] = SAM[i, "GOV"] - (Xv00[i] - SAM[i, "INV"])
# Computing the direct tax revenue that balances the gov. budget
Td00[...] = Sum(i, Xg00[i]) - Sum(i, Tz00[i] + Tm00[i])
# Computing the household sav. that balances the household budget
Sp00[...] = Sum(h, FF00[h]) - (Sum(i, Xp00[i]) + Td00)
III00[...] = Sum(i, Xv00[i])
II00[j] = (Sp00 + Sf00) * F00["CAP", j] / Sum(i, F00["CAP", i])
KK00[j] = F00["CAP", j] / ror
# ===============================================================
# Computing the BAU path
# ===============================================================
Y0[j, t] = Y00[j] * (1 + pop) ** (Ord(t) - 1)
F0[h, j, t] = F00[h, j] * (1 + pop) ** (Ord(t) - 1)
X0[i, j, t] = X00[i, j] * (1 + pop) ** (Ord(t) - 1)
Z0[j, t] = Z00[j] * (1 + pop) ** (Ord(t) - 1)
Xp0[i, t] = Xp00[i] * (1 + pop) ** (Ord(t) - 1)
Xv0[i, t] = Xv00[i] * (1 + pop) ** (Ord(t) - 1)
E0[i, t] = E00[i] * (1 + pop) ** (Ord(t) - 1)
M0[i, t] = M00[i] * (1 + pop) ** (Ord(t) - 1)
Q0[i, t] = Q00[i] * (1 + pop) ** (Ord(t) - 1)
D0[i, t] = D00[i] * (1 + pop) ** (Ord(t) - 1)
FF0[h, t] = FF00[h] * (1 + pop) ** (Ord(t) - 1)
III0[t] = III00 * (1 + pop) ** (Ord(t) - 1)
II0[j, t] = II00[j] * (1 + pop) ** (Ord(t) - 1)
KK0[j, t] = KK00[j] * (1 + pop) ** (Ord(t) - 1)
CC0[t] = CC00 * (1 + pop) ** (Ord(t) - 1)
Sp0[t] = Sp00 * (1 + pop) ** (Ord(t) - 1)
Td0[t] = Td00 * (1 + pop) ** (Ord(t) - 1)
Tz0[j, t] = Tz00[j] * (1 + pop) ** (Ord(t) - 1)
Tm0[i, t] = Tm00[i] * (1 + pop) ** (Ord(t) - 1)
pf0[h, j, t] = 1
py0[j, t] = 1
pz0[j, t] = 1
pq0[i, t] = 1
pe0[i, t] = 1
pm0[i, t] = 1
pd0[i, t] = 1
pk0[t] = 1
epsilon0[t] = 1
PRICE0[t] = 1
# Setting exogenous variables
Xg0[i, t] = Xg00[i] * (1 + pop) ** (Ord(t) - 1)
Sf0[t] = Sf00 * (1 + pop) ** (Ord(t) - 1)
pWe[i] = 1
pWm[i] = 1
tauz[i] = tauz00[i]
taum[i] = taum00[i]
# ===============================================================
# Calibration ---------------------------------------------------
# ===============================================================
# Parameters
sigma = Parameter(
m, name="sigma", domain=i, description="elasticity of substitution"
)
psi = Parameter(
m, name="psi", domain=i, description="elasticity of transformation"
)
eta = Parameter(
m,
name="eta",
domain=i,
description="substitution elasticity parameter",
)
phi = Parameter(
m,
name="phi",
domain=i,
description="transformation elasticity parameter",
)
sigma[i] = 2
psi[i] = 2
eta[i] = (sigma[i] - 1) / sigma[i]
phi[i] = (psi[i] + 1) / psi[i]
# Parameters
alpha = Parameter(
m,
name="alpha",
domain=i,
description="share par. in composite cons. func.",
)
a = Parameter(
m, name="a", description="scale par. in composite cons. func."
)
beta = Parameter(
m,
name="beta",
domain=[h, j],
description="share par. in production func.",
)
b = Parameter(
m, name="b", domain=j, description="scale par. in production func."
)
ax = Parameter(
m,
name="ax",
domain=[i, j],
description="intermediate input requirement coeff.",
)
ay = Parameter(
m, name="ay", domain=j, description="composite fact. input req. coeff."
)
lamda = Parameter(
m, name="lamda", domain=i, description="investment demand share"
)
iota = Parameter(
m, name="iota", description="scale par. in comp. inv. prod. func."
)
deltam = Parameter(
m, name="deltam", domain=i, description="share par. in Armington func."
)
deltad = Parameter(
m, name="deltad", domain=i, description="share par. in Armington func."
)
gamma = Parameter(
m, name="gamma", domain=i, description="scale par. in Armington func."
)
xid = Parameter(
m,
name="xid",
domain=i,
description="share par. in transformation func.",
)
xie = Parameter(
m,
name="xie",
domain=i,
description="share par. in transformation func.",
)
theta = Parameter(
m,
name="theta",
domain=i,
description="scale par. in transformation func.",
)
ssp = Parameter(m, name="ssp", description="propensity to save")
alpha[i] = Xp00[i] / Sum(j, Xp00[j])
a[...] = CC00 / Product(j, Xp00[j] ** alpha[j])
beta[h, j] = F00[h, j] / Sum(k, F00[k, j])
b[j] = Y00[j] / Product(h, F00[h, j] ** beta[h, j])
ax[i, j] = X00[i, j] / Z00[j]
ay[j] = Y00[j] / Z00[j]
lamda[i] = Xv00[i] / Sum(j, Xv00[j])
iota[...] = III00 / Product(i, Xv00[i] ** lamda[i])
deltam[i] = (
(1 + taum00[i])
* M00[i] ** (1 - eta[i])
/ ((1 + taum00[i]) * M00[i] ** (1 - eta[i]) + D00[i] ** (1 - eta[i]))
)
deltad[i] = D00[i] ** (1 - eta[i]) / (
(1 + taum00[i]) * M00[i] ** (1 - eta[i]) + D00[i] ** (1 - eta[i])
)
gamma[i] = Q00[i] / (
deltam[i] * M00[i] ** eta[i] + deltad[i] * D00[i] ** eta[i]
) ** (1 / eta[i])
xie[i] = E00[i] ** (1 - phi[i]) / (
E00[i] ** (1 - phi[i]) + D00[i] ** (1 - phi[i])
)
xid[i] = D00[i] ** (1 - phi[i]) / (
E00[i] ** (1 - phi[i]) + D00[i] ** (1 - phi[i])
)
theta[i] = Z00[i] / (
xie[i] * E00[i] ** phi[i] + xid[i] * D00[i] ** phi[i]
) ** (1 / phi[i])
ssp[...] = Sp00 / (Sum([h, j], F00[h, j]) - Td00)
# ===============================================================
# Defining model system -----------------------------------------
# ===============================================================
# Variables
Y = Variable(
m, name="Y", type="free", domain=j, description="composite factor"
)
F = Variable(
m, name="F", type="free", domain=[h, j], description="factor input"
)
X = Variable(
m,
name="X",
type="free",
domain=[i, j],
description="intermediate input",
)
Z = Variable(
m, name="Z", type="free", domain=j, description="gross domestic output"
)
Xp = Variable(
m,
name="Xp",
type="free",
domain=i,
description="household consumption",
)
Xg = Variable(
m,
name="Xg",
type="free",
domain=i,
description="government consumption",
)
Xv = Variable(
m, name="Xv", type="free", domain=i, description="investment demand"
)
E = Variable(m, name="E", type="free", domain=i, description="exports")
M = Variable(m, name="M", type="free", domain=i, description="imports")
Q = Variable(
m,
name="Q",
type="free",
domain=i,
description="Armington's composite good",
)
D = Variable(
m, name="D", type="free", domain=i, description="domestic good"
)
FF = Variable(
m, name="FF", type="free", domain=h, description="factor endowments"
)
pf = Variable(
m, name="pf", type="free", domain=[h, j], description="factor price"
)
py = Variable(
m,
name="py",
type="free",
domain=j,
description="composite factor price",
)
pz = Variable(
m,
name="pz",
type="free",
domain=j,
description="supply price of gross domestic output",
)
pq = Variable(
m,
name="pq",
type="free",
domain=i,
description="Armington's composite good price",
)
pe = Variable(
m,
name="pe",
type="free",
domain=i,
description="export price in local currency",
)
pm = Variable(
m,
name="pm",
type="free",
domain=i,
description="import price in local currency",
)
pd = Variable(
m, name="pd", type="free", domain=i, description="domestic good price"
)
pk = Variable(
m,
name="pk",
type="free",
description="composite investment goods price",
)
epsilon = Variable(
m, name="epsilon", type="free", description="exchange rate"
)
Sp = Variable(m, name="Sp", type="free", description="private savings")
Sf = Variable(m, name="Sf", type="free", description="foreign savings")
Td = Variable(m, name="Td", type="free", description="direct tax")
Tz = Variable(
m, name="Tz", type="free", domain=j, description="production tax"
)
Tm = Variable(
m, name="Tm", type="free", domain=i, description="import tariff"
)
KK = Variable(
m, name="KK", type="free", domain=j, description="capital stock"
)
II = Variable(
m, name="II", type="free", domain=j, description="sectoral investment"
)
III = Variable(
m, name="III", type="free", description="composite investment"
)
PRICE = Variable(
m, name="PRICE", type="free", description="numeraire price"
)
# Equations
eqpy = Equation(
m, name="eqpy", domain=j, description="composite factor prod. func."
)
eqF = Equation(
m, name="eqF", domain=[h, j], description="factor demand function"
)
eqX = Equation(
m,
name="eqX",
domain=[i, j],
description="intermediate demand function",
)
eqY = Equation(
m, name="eqY", domain=j, description="composite factor demand function"
)
eqpzs = Equation(
m, name="eqpzs", domain=j, description="unit cost function"
)
eqTd = Equation(m, name="eqTd", description="direct tax revenue function")
eqTz = Equation(
m, name="eqTz", domain=j, description="production tax revenue function"
)
eqTm = Equation(
m, name="eqTm", domain=i, description="import tariff revenue function"
)
eqXv = Equation(
m, name="eqXv", domain=i, description="investment demand function"
)
eqSp = Equation(m, name="eqSp", description="private saving function")
eqXp = Equation(
m, name="eqXp", domain=i, description="household demand function"
)
eqpe = Equation(
m, name="eqpe", domain=i, description="world export price equation"
)
eqpm = Equation(
m, name="eqpm", domain=i, description="world import price equation"
)
eqepsilon = Equation(
m, name="eqepsilon", description="balance of payments"
)
eqpqs = Equation(
m, name="eqpqs", domain=i, description="Armington function"
)
eqM = Equation(
m, name="eqM", domain=i, description="import demand function"
)
eqD = Equation(
m, name="eqD", domain=i, description="domestic good demand function"
)
eqpzd = Equation(
m, name="eqpzd", domain=i, description="transformation function"
)
eqE = Equation(
m, name="eqE", domain=i, description="export supply function"
)
eqDs = Equation(
m, name="eqDs", domain=i, description="domestic good supply function"
)
eqpqd = Equation(
m,
name="eqpqd",
domain=i,
description="market clearing cond. for comp. good",
)
eqpf1 = Equation(
m,
name="eqpf1",
domain=[h_mob],
description="mobile factor market clearing cond.",
)
eqpf2 = Equation(
m,
name="eqpf2",
domain=[h_mob, i, j],
description="mobile factor market clearing cond.",
)
eqpf3 = Equation(
m,
name="eqpf3",
domain=j,
description="immobile factor market clearing cond.",
)
eqpk = Equation(
m, name="eqpk", description="composite inv. goods mar. clear. cond."
)
eqIII = Equation(
m, name="eqIII", description="composite inv. goods production func."
)
eqII = Equation(
m,
name="eqII",
domain=j,
description="evolution of target capital stocks",
)
eqPRICE = Equation(m, name="eqPRICE", description="numeraire price")
# ===============================================================
# Model equations
# ===============================================================
# [domestic production] -
# composite factor production func. (Cobb-Douglas)
eqpy[j] = Y[j] == b[j] * Product(h, F[h, j] ** beta[h, j])
# factor demand function (Cobb-Douglas)
eqF[h, j] = F[h, j] == beta[h, j] * py[j] * Y[j] / pf[h, j]
# intermediate input demand function (Leontief)
eqX[i, j] = X[i, j] == ax[i, j] * Z[j]
# composite factor demand function (Leontief)
eqY[j] = Y[j] == ay[j] * Z[j]
# unit price of gross output (Leontief)
eqpzs[j] = pz[j] == ay[j] * py[j] + Sum(i, ax[i, j] * pq[i])
# [government behavior] -
# lump Sum direct tax revenue
eqTd[...] = Td == Sum(i, pq[i] * Xg[i]) - Sum(i, Tm[i] + Tz[i])
# production tax revenue
eqTz[j] = Tz[j] == tauz[j] * pz[j] * Z[j]
# import tariff revenue
eqTm[i] = Tm[i] == taum[i] * pm[i] * M[i]
# [investment behavior] -
# composite investment production function
eqXv[i] = Xv[i] == lamda[i] * pk * Sum(j, II[j]) / pq[i]
# [savings] ----------
# savings function
eqSp[...] = Sp == ssp * (Sum([h, j], pf[h, j] * F[h, j]) - Td)
# [household consumption] -- (Cobb-Douglas)
eqXp[i] = (
Xp[i] == alpha[i] * (Sum([h, j], pf[h, j] * F[h, j]) - Sp - Td) / pq[i]
)
# [international trade] --
eqpe[i] = pe[i] == epsilon * pWe[i]
eqpm[i] = pm[i] == epsilon * pWm[i]
# BOP constraint
eqepsilon[...] = Sum(i, pWe[i] * E[i]) + Sf == Sum(i, pWm[i] * M[i])
# [Armington function] --
# Armington's composite good production function (CES)
eqpqs[i] = Q[i] == gamma[i] * (
deltam[i] * M[i] ** eta[i] + deltad[i] * D[i] ** eta[i]
) ** (1 / eta[i])
# import demand function (CES)
eqM[i] = (
M[i]
== (gamma[i] ** eta[i] * deltam[i] * pq[i] / ((1 + taum[i]) * pm[i]))
** (1 / (1 - eta[i]))
* Q[i]
)
# domestic good demand function (CES)
eqD[i] = (
D[i]
== (gamma[i] ** eta[i] * deltad[i] * pq[i] / pd[i])
** (1 / (1 - eta[i]))
* Q[i]
)
# [transformation function] --
# gross domestic output disaggregation function (CET)
eqpzd[i] = Z[i] == theta[i] * (
xie[i] * E[i] ** phi[i] + xid[i] * D[i] ** phi[i]
) ** (1 / phi[i])
# export supply function (CET)
eqE[i] = (
E[i]
== (theta[i] ** phi[i] * xie[i] * (1 + tauz[i]) * pz[i] / pe[i])
** (1 / (1 - phi[i]))
* Z[i]
)
# domestic good supply function (CET)
eqDs[i] = (
D[i]
== (theta[i] ** phi[i] * xid[i] * (1 + tauz[i]) * pz[i] / pd[i])
** (1 / (1 - phi[i]))
* Z[i]
)
# [market clearing condition]
# Arminton's composite good market
eqpqd[i] = Q[i] == Xp[i] + Xg[i] + Xv[i] + Sum(j, X[i, j])
# labor market: quantity
eqpf1[h_mob] = Sum(j, F[h_mob, j]) == FF[h_mob]
# labor market: price
eqpf2[h_mob, i, j] = pf[h_mob, j] == pf[h_mob, i]
# capital market
eqpf3[j] = F["CAP", j] == ror * KK[j]
# investment goods market
eqpk[...] = Sum(j, II[j]) == III
# [dynamic equations]
# composite investment good market clearing condition
eqIII[...] = iota * Product(i, Xv[i] ** lamda[i]) == III
# sectoral investment allocation
eqII[j] = pk * II[j] == pf["CAP", j] ** zeta * F["CAP", j] / Sum(
i, pf["CAP", i] ** zeta * F["CAP", i]
) * (Sp + epsilon * Sf)
# felicity function
CC = a * Product(i, Xp[i] ** alpha[i])
# Price level [numeraire]
eqPRICE[...] = Sum(j, pq[j] * Q00[j] / Sum(i, Q00[i])) == PRICE
# ===============================================================
# Initializing variables ----------------------------------------
# ===============================================================
Y.l[j] = Y00[j]
F.l[h, j] = F00[h, j]
X.l[i, j] = X00[i, j]
Z.l[j] = Z00[j]
Xp.l[i] = Xp00[i]
Xv.l[i] = Xv00[i]
E.l[i] = E00[i]
M.l[i] = M00[i]
Q.l[i] = Q00[i]
D.l[i] = D00[i]
pf.l[h, j] = 1
py.l[j] = 1
pz.l[j] = 1
pq.l[i] = 1
pe.l[i] = 1
pm.l[i] = 1
pd.l[i] = 1
pk.l[...] = 1
epsilon.l[...] = 1
Sp.l[...] = Sp00
Td.l[...] = Td00
Tz.l[j] = Tz00[j]
Tm.l[i] = Tm00[i]
FF.l[h] = FF00[h]
III.l[...] = III00
II.l[j] = II00[j]
# ---------------------------------------------------------------
# Numeraire
PRICE.fx[...] = 1
# Initial factor endowments and exogenous variables
FF.fx[h_mob] = FF00[h_mob]
KK.fx[j] = KK00[j]
Xg.fx[i] = Xg00[i]
Sf.fx[...] = Sf00
# ===============================================================
# Defining and solving the model --------------------------------
# ===============================================================
dyncge = Model(
m,
name="dyncge",
equations=m.getEquations(),
problem=Problem.NLP,
sense=Sense.MAX,
objective=CC,
)
dyncge.solve()
# ===============================================================
# Simulation Runs: Abolition of Import Tariffs
# ===============================================================
# Scenario:
taum[i] = taum00[i] * 0
for iteration, _ in t.records.itertuples(index=False):
dyncge.solve()
# storing results -------------------------
Y1[j, iteration] = Y.l[j]
F1[h, j, iteration] = F.l[h, j]
X1[i, j, iteration] = X.l[i, j]
Z1[j, iteration] = Z.l[j]
Xp1[i, iteration] = Xp.l[i]
Xv1[i, iteration] = Xv.l[i]
E1[i, iteration] = E.l[i]
M1[i, iteration] = M.l[i]
Q1[i, iteration] = Q.l[i]
D1[i, iteration] = D.l[i]
Sp1[iteration] = Sp.l
Td1[iteration] = Td.l
Tz1[j, iteration] = Tz.l[j]
Tm1[i, iteration] = Tm.l[i]
FF1[h, iteration] = FF.l[h]
II1[j, iteration] = II.l[j]
III1[iteration] = III.l
KK1[j, iteration] = KK.l[j]
CC1[iteration] = dyncge.objective_value
tauz1[i, iteration] = tauz[i]
taum1[i, iteration] = taum[i]
pf1[h, j, iteration] = pf.l[h, j]
py1[j, iteration] = py.l[j]
pz1[j, iteration] = pz.l[j]
pd1[j, iteration] = pd.l[j]
pe1[j, iteration] = pe.l[j]
pm1[j, iteration] = pm.l[j]
pq1[j, iteration] = pq.l[j]
pk1[iteration] = pk.l
epsilon1[iteration] = epsilon.l
PRICE1[iteration] = PRICE.l
# updating the state variables --------------
FF.fx[h_mob] = FF.l[h_mob] * (1 + pop)
KK.fx[j] = (1 - dep) * KK.l[j] + II.l[j]
if int(iteration) < 30:
Xg.fx[i] = Xg0[i, str(int(iteration) + 1)]
Sf.fx[...] = Sf0[str(int(iteration) + 1)]
# ===============================================================
# Aftermath Computation
# ===============================================================
# Display of changes --------------------------------------------
import math
assert math.isclose(dyncge.objective_value, 539570.5027, rel_tol=0.001)
# Parameters
# changes
dY = Parameter(
m,
name="dY",
domain=[j, t],
description="change of composite factor [%]",
)
dF = Parameter(
m,
name="dF",
domain=[h, j, t],
description="change of factor input [%]",
)
dX = Parameter(
m,
name="dX",
domain=[i, j, t],
description="change of intermediate input [%]",
)
dZ = Parameter(
m,
name="dZ",
domain=[j, t],
description="change of gross output [%]",
)
dXp = Parameter(
m,
name="dXp",
domain=[i, t],
description="change of household consumption [%]",
)
dXv = Parameter(
m,
name="dXv",
domain=[i, t],
description="change of investment demand [%]",
)
dE = Parameter(
m,
name="dE",
domain=[i, t],
description="change of exports [%]",
)
dM = Parameter(
m,
name="dM",
domain=[i, t],
description="change of imports [%]",
)
dQ = Parameter(
m,
name="dQ",
domain=[i, t],
description="change of Armington's composite good [%]",
)
dD = Parameter(
m,
name="dD",
domain=[i, t],
description="change of domestic good [%]",
)
dSp = Parameter(
m,
name="dSp",
domain=t,
description="change of private saving [%]",
)
dTd = Parameter(
m,
name="dTd",
domain=t,
description="change of direct tax [%]",
)
dTz = Parameter(
m,
name="dTz",
domain=[j, t],
description="change of production tax [%]",
)
dTm = Parameter(
m,
name="dTm",
domain=[i, t],
description="change of import tariff [%]",
)
dFF = Parameter(
m,
name="dFF",
domain=[h, t],
description="change of initial sectoral factor uses [%]",
)
dKK = Parameter(
m,
name="dKK",
domain=[j, t],
description="change of sectoral capital stock [%]",
)
dII = Parameter(
m,
name="dII",
domain=[j, t],
description="change of sectoral investment [%]",
)
dIII = Parameter(
m,
name="dIII",
domain=t,
description="change of composite investment [%]",
)
dCC = Parameter(
m,
name="dCC",
domain=t,
description="change of utility [%]",
)
dpz = Parameter(
m,
name="dpz",
domain=[j, t],
description="change of gross output price [%]",
)
dpd = Parameter(
m,
name="dpd",
domain=[j, t],
description="change of domestic good price [%]",
)
dpm = Parameter(
m,
name="dpm",
domain=[j, t],
description="change of import price [%]",
)
dpe = Parameter(
m,
name="dpe",
domain=[j, t],
description="change of export price [%]",
)
dpq = Parameter(
m,
name="dpq",
domain=[j, t],
description="change of Armington's comp. good price [%]",
)
dpf = Parameter(
m,
name="dpf",
domain=[h, j, t],
description="change of factor price [%]",
)
dpy = Parameter(
m,
name="dpy",
domain=[j, t],
description="change of composite factor price [%]",
)
depsilon = Parameter(
m,
name="depsilon",
domain=t,
description="change of foreign exchange rate [%]",
)
dpk = Parameter(
m,
name="dpk",
domain=t,
description="change of capital good price [%]",
)
# BAU growth rate
gY0 = Parameter(
m,
name="gY0",
domain=[j, t],
description="growth of composite factor [%]",
)
gF0 = Parameter(
m,
name="gF0",
domain=[h, j, t],
description="growth of factor input [%]",
)
gX0 = Parameter(
m,
name="gX0",
domain=[i, j, t],
description="growth of intermediate input [%]",
)
gZ0 = Parameter(
m,
name="gZ0",
domain=[j, t],
description="growth of gross output [%]",
)
gXp0 = Parameter(
m,
name="gXp0",
domain=[i, t],
description="growth of household consumption [%]",
)
gXv0 = Parameter(
m,
name="gXv0",
domain=[i, t],
description="growth of investment demand [%]",
)
gE0 = Parameter(
m,
name="gE0",
domain=[i, t],
description="growth of exports [%]",
)
gM0 = Parameter(
m,
name="gM0",
domain=[i, t],
description="growth of imports [%]",
)
gQ0 = Parameter(
m,
name="gQ0",
domain=[i, t],
description="growth of Armington's composite good [%]",
)
gD0 = Parameter(
m,
name="gD0",
domain=[i, t],
description="growth of domestic good [%]",
)
gSp0 = Parameter(
m,
name="gSp0",
domain=t,
description="growth of private saving [%]",
)
gTd0 = Parameter(
m,
name="gTd0",
domain=t,
description="growth of direct tax [%]",
)
gTz0 = Parameter(
m,
name="gTz0",
domain=[j, t],
description="growth of production tax [%]",
)
gTm0 = Parameter(
m,
name="gTm0",
domain=[i, t],
description="growth of import tariff [%]",
)
gFF0 = Parameter(
m,
name="gFF0",
domain=[h, t],
description="growth of initial sectoral factor uses [%]",
)
gKK0 = Parameter(
m,
name="gKK0",
domain=[j, t],
description="growth of sectoral capital stock [%]",
)
gII0 = Parameter(
m,
name="gII0",
domain=[j, t],
description="growth of sectoral investment [%]",
)
gIII0 = Parameter(
m,
name="gIII0",
domain=t,
description="growth of composite investment [%]",
)
gCC0 = Parameter(
m,
name="gCC0",
domain=t,
description="growth of growth rate of CC [%]",
)
# C/F growth rate
gY1 = Parameter(
m,
name="gY1",
domain=[j, t],
description="growth of composite factor [%]",
)
gF1 = Parameter(
m,
name="gF1",
domain=[h, j, t],
description="growth of factor input [%]",
)
gX1 = Parameter(
m,
name="gX1",
domain=[i, j, t],
description="growth of intermediate input [%]",
)
gZ1 = Parameter(
m,
name="gZ1",
domain=[j, t],
description="growth of gross output [%]",
)
gXp1 = Parameter(
m,
name="gXp1",
domain=[i, t],
description="growth of household consumption [%]",
)
gXv1 = Parameter(
m,
name="gXv1",
domain=[i, t],
description="growth of investment demand [%]",
)
gE1 = Parameter(
m,
name="gE1",
domain=[i, t],
description="growth of exports [%]",
)
gM1 = Parameter(
m,
name="gM1",
domain=[i, t],
description="growth of imports [%]",
)
gQ1 = Parameter(
m,
name="gQ1",
domain=[i, t],
description="growth of Armington's composite good [%]",
)
gD1 = Parameter(
m,
name="gD1",
domain=[i, t],
description="growth of domestic good [%]",
)
gSp1 = Parameter(
m,
name="gSp1",
domain=t,
description="growth of private saving [%]",
)
gTd1 = Parameter(
m,
name="gTd1",
domain=t,
description="growth of direct tax [%]",
)
gTz1 = Parameter(
m,
name="gTz1",
domain=[j, t],
description="growth of production tax [%]",
)
gTm1 = Parameter(
m,
name="gTm1",
domain=[i, t],
description="growth of import tariff [%]",
)
gFF1 = Parameter(
m,
name="gFF1",
domain=[h, t],
description="growth of initial sectoral factor uses [%]",
)
gKK1 = Parameter(
m,
name="gKK1",
domain=[j, t],
description="growth of sectoral capital stock [%]",
)
gII1 = Parameter(
m,
name="gII1",
domain=[j, t],
description="growth of sectoral investment [%]",
)
gIII1 = Parameter(
m,
name="gIII1",
domain=t,
description="growth of composite investment [%]",
)
gCC1 = Parameter(
m,
name="gCC1",
domain=t,
description="growth of growth rate of CC [%]",
)
# welfare
EV = Parameter(
m, name="EV", domain=t, description="equivalent variations [current]"
)
EV_TTL = Parameter(
m, name="EV_TTL", description="total EV [discounted Sum]"
)
dY[j, t].where[Y0[j, t]] = (Y1[j, t] / Y0[j, t] - 1) * 100
dF[h, j, t].where[F0[h, j, t]] = (F1[h, j, t] / F0[h, j, t] - 1) * 100
dX[i, j, t].where[X0[i, j, t]] = (X1[i, j, t] / X0[i, j, t] - 1) * 100
dZ[j, t].where[Z0[j, t]] = (Z1[j, t] / Z0[j, t] - 1) * 100
dXp[i, t].where[Xp0[i, t]] = (Xp1[i, t] / Xp0[i, t] - 1) * 100
dXv[i, t].where[Xv0[i, t]] = (Xv1[i, t] / Xv0[i, t] - 1) * 100
dE[i, t].where[E0[i, t]] = (E1[i, t] / E0[i, t] - 1) * 100
dM[i, t].where[M0[i, t]] = (M1[i, t] / M0[i, t] - 1) * 100
dQ[i, t].where[Q0[i, t]] = (Q1[i, t] / Q0[i, t] - 1) * 100
dD[i, t].where[D0[i, t]] = (D1[i, t] / D0[i, t] - 1) * 100
dSp[t].where[Sp0[t]] = (Sp1[t] / Sp0[t] - 1) * 100
dTd[t].where[Td0[t]] = (Td1[t] / Td0[t] - 1) * 100
dTz[j, t].where[Tz0[j, t]] = (Tz1[j, t] / Tz0[j, t] - 1) * 100
dTm[i, t].where[Tm0[i, t]] = (Tm1[i, t] / Tm0[i, t] - 1) * 100
dFF[h, t].where[FF0[h, t]] = (FF1[h, t] / FF0[h, t] - 1) * 100
dII[j, t].where[II0[j, t]] = (II1[j, t] / II0[j, t] - 1) * 100
dIII[t].where[III0[t]] = (III1[t] / III0[t] - 1) * 100
dKK[j, t].where[KK0[j, t]] = (KK1[j, t] / KK0[j, t] - 1) * 100
dCC[t].where[CC0[t]] = (CC1[t] / CC0[t] - 1) * 100
dpz[j, t].where[pz0[j, t]] = (pz1[j, t] / pz0[j, t] - 1) * 100
dpd[j, t].where[pd0[j, t]] = (pd1[j, t] / pd0[j, t] - 1) * 100
dpm[j, t].where[pm0[j, t]] = (pm1[j, t] / pm0[j, t] - 1) * 100
dpe[j, t].where[pe0[j, t]] = (pe1[j, t] / pe0[j, t] - 1) * 100
dpq[j, t].where[pq0[j, t]] = (pq1[j, t] / pq0[j, t] - 1) * 100
dpf[h, j, t].where[pf0[h, j, t]] = (pf1[h, j, t] / pf0[h, j, t] - 1) * 100
dpy[j, t].where[py0[j, t]] = (py1[j, t] / py0[j, t] - 1) * 100
depsilon[t].where[epsilon0[t]] = (epsilon1[t] / epsilon0[t] - 1) * 100
dpk[t].where[pk0[t]] = (pk1[t] / pk0[t] - 1) * 100
gY0[j, t.lead(1)].where[Y0[j, t]] = (Y0[j, t.lead(1)] / Y0[j, t] - 1) * 100
gF0[h, j, t.lead(1)].where[F0[h, j, t]] = (
F0[h, j, t.lead(1)] / F0[h, j, t] - 1
) * 100
gX0[i, j, t.lead(1)].where[X0[i, j, t]] = (
X0[i, j, t.lead(1)] / X0[i, j, t] - 1
) * 100
gZ0[j, t.lead(1)].where[Z0[j, t]] = (Z0[j, t.lead(1)] / Z0[j, t] - 1) * 100
gXp0[i, t.lead(1)].where[Xp0[i, t]] = (
Xp0[i, t.lead(1)] / Xp0[i, t] - 1
) * 100
gXv0[i, t.lead(1)].where[Xv0[i, t]] = (
Xv0[i, t.lead(1)] / Xv0[i, t] - 1
) * 100
gE0[i, t.lead(1)].where[E0[i, t]] = (E0[i, t.lead(1)] / E0[i, t] - 1) * 100
gM0[i, t.lead(1)].where[M0[i, t]] = (M0[i, t.lead(1)] / M0[i, t] - 1) * 100
gQ0[i, t.lead(1)].where[Q0[i, t]] = (Q0[i, t.lead(1)] / Q0[i, t] - 1) * 100
gD0[i, t.lead(1)].where[D0[i, t]] = (D0[i, t.lead(1)] / D0[i, t] - 1) * 100
gSp0[t.lead(1)].where[Sp0[t]] = (Sp0[t.lead(1)] / Sp0[t] - 1) * 100
gTd0[t.lead(1)].where[Td0[t]] = (Td0[t.lead(1)] / Td0[t] - 1) * 100
gTz0[j, t.lead(1)].where[Tz0[j, t]] = (
Tz0[j, t.lead(1)] / Tz0[j, t] - 1
) * 100
gTm0[i, t.lead(1)].where[Tm0[i, t]] = (
Tm0[i, t.lead(1)] / Tm0[i, t] - 1
) * 100
gFF0[h, t.lead(1)].where[FF0[h, t]] = (
FF0[h, t.lead(1)] / FF0[h, t] - 1
) * 100
gII0[j, t.lead(1)].where[II0[j, t]] = (
II0[j, t.lead(1)] / II0[j, t] - 1
) * 100
gIII0[t.lead(1)].where[III0[t]] = (III0[t.lead(1)] / III0[t] - 1) * 100
gKK0[j, t.lead(1)].where[KK0[j, t]] = (
KK0[j, t.lead(1)] / KK0[j, t] - 1
) * 100
gCC0[t.lead(1)].where[CC0[t]] = (CC0[t.lead(1)] / CC0[t] - 1) * 100
gY1[j, t.lead(1)].where[Y1[j, t]] = (Y1[j, t.lead(1)] / Y1[j, t] - 1) * 100
gF1[h, j, t.lead(1)].where[F1[h, j, t]] = (
F1[h, j, t.lead(1)] / F1[h, j, t] - 1
) * 100
gX1[i, j, t.lead(1)].where[X1[i, j, t]] = (
X1[i, j, t.lead(1)] / X1[i, j, t] - 1
) * 100
gZ1[j, t.lead(1)].where[Z1[j, t]] = (Z1[j, t.lead(1)] / Z1[j, t] - 1) * 100
gXp1[i, t.lead(1)].where[Xp1[i, t]] = (
Xp1[i, t.lead(1)] / Xp1[i, t] - 1
) * 100
gXv1[i, t.lead(1)].where[Xv1[i, t]] = (
Xv1[i, t.lead(1)] / Xv1[i, t] - 1
) * 100
gE1[i, t.lead(1)].where[E1[i, t]] = (E1[i, t.lead(1)] / E1[i, t] - 1) * 100
gM1[i, t.lead(1)].where[M1[i, t]] = (M1[i, t.lead(1)] / M1[i, t] - 1) * 100
gQ1[i, t.lead(1)].where[Q1[i, t]] = (Q1[i, t.lead(1)] / Q1[i, t] - 1) * 100
gD1[i, t.lead(1)].where[D1[i, t]] = (D1[i, t.lead(1)] / D1[i, t] - 1) * 100
gSp1[t.lead(1)].where[Sp1[t]] = (Sp1[t.lead(1)] / Sp1[t] - 1) * 100
gTd1[t.lead(1)].where[Td1[t]] = (Td1[t.lead(1)] / Td1[t] - 1) * 100
gTz1[j, t.lead(1)].where[Tz1[j, t]] = (
Tz1[j, t.lead(1)] / Tz1[j, t] - 1
) * 100
gTm1[i, t.lead(1)].where[Tm1[i, t]] = (
Tm1[i, t.lead(1)] / Tm1[i, t] - 1
) * 100
gFF1[h, t.lead(1)].where[FF1[h, t]] = (
FF1[h, t.lead(1)] / FF1[h, t] - 1
) * 100
gII1[j, t.lead(1)].where[II1[j, t]] = (
II1[j, t.lead(1)] / II1[j, t] - 1
) * 100
gIII1[t.lead(1)].where[III1[t]] = (III1[t.lead(1)] / III1[t] - 1) * 100
gKK1[j, t.lead(1)].where[KK1[j, t]] = (
KK1[j, t.lead(1)] / KK1[j, t] - 1
) * 100
gCC1[t.lead(1)].where[CC1[t]] = (CC1[t.lead(1)] / CC1[t] - 1) * 100
# Welfare measure: Hicksian equivalent variations ---------------
EV[t] = (CC1[t] - CC0[t]) / a / Product(i, (alpha[i] / 1) ** alpha[i])
EV_TTL[...] = Sum(t, EV[t] / (1 + ror) ** (Ord(t) - 1))
print("EV_TTL: ", round(EV_TTL.records.value[0], 3))
if __name__ == "__main__":
main()
