Introducing conductance as our production capacity and free energy as the driving force, we can create a matrix of velocities, which at first seems unnecessary but becomes incredibly useful later.

# Introducing conductance x free energy factor
N = [100, 0]
t =
0.0001
N_hst = [N]
sig_x_A = sxA = [[
0, 2], [1, 0]] # Defining our sigma_jk * A_jk matrix
for i in range(0, 25000):
    v = [sxA[
0][1]*N[1] - sxA[1][0]*N[0], sxA[1][0]*N[0] - sxA[0][1]*N[1]]
   
# Exactly the same as before but by pulling the values from sxA
   
N = [N[0] + v[0] * t, N[1] + v[1] * t]
    N_hst.append(N)
plt.plot(
range(0, len(N_hst)), [i[0] for i in N_hst])
plt.plot(
range(0, len(N_hst)), [i[1] for i in N_hst])
plt.show()

Now we have a ‘sigma A’, denoted sxA, matrix which, if we take ‘k_B T=1’, is the ‘v_jk’ (eq.5) matrix. Now we can observe that the matrix laid out in the code above is equivalent to the ratio we introduced by hand earlier, due to the identical graphs. We should easily be able to introduce more commodities with this improved model. It important to realise the inclusion of N in the velocities at this point. The free energy is per ‘quantity’ of the commodity so to find the total driving force behind a given commodity requires scaling by its current worth/amount.

# Introducing a third commodity
N = [50, 10, 40] # Added a third starting value
t = 0.0001
N_hst = [N]
sig_x_A = sxA = [[
0, 1, 1], # Now a 3x3 matrix is required to describe
   
[1, 0, 1], # all possible relationships (self-self are 0)
   
[1, 1, 0]]
for i in range(0, 25000):
    v = [sxA[
0][1] * N[1] + sxA[0][2] * N[2] - sxA[1][0] * N[0] - sxA[2][0] * N[0],
   
sxA[1][0] * N[0] + sxA[1][2] * N[2] - sxA[0][1] * N[1] - sxA[2][1] * N[1],
   
sxA[2][0] * N[0] + sxA[2][1] * N[1] - sxA[0][2] * N[2] - sxA[1][2] * N[2]]
   
# Two new terms per existing commodity for each new one
   
N = [N[0] + v[0] * t,
   
N[1] + v[1] * t,
   
N[2] + v[2] * t]
   
# A new line to update each new commodity
   
N_hst.append(N)
plt.plot(
range(0, len(N_hst)), [i[0] for i in N_hst])
plt.plot(
range(0, len(N_hst)), [i[1] for i in N_hst])
plt.plot(
range(0, len(N_hst)), [i[2] for i in N_hst])
plt.show()

The ability to track three commodities is clearly successful, they all tend to equal worth in the figure when they are valued equally, but it requires a large manual input to track a new commodity. The velocities, v, and updating N, should be done based on input.

# Identify and implement emergent pattern in velocity updating
N = [50, 10, 40]
t =
0.0001
N_hst = [N]
sig_x_A = sxA = [[
0, 1, 1],
   
[1, 0, 1],
   
[1, 1, 0]]
''' v indexed jk, -> v_jk:
notation sxA[0][1]*N[1] -> 011
v_j = (j00 + j11 + j22 ... + jkk) - (0jj + 1jj + 2jj ... + kjj)'''
def v_j(j): # A function that determines the appropriate v
    v_ =
0
   
for k in range(0, len(N)): # Repeats to cover all combinations
        v_ += sxA[j][k] * N[k]
# Annila & Salthe velocity
        v_ += -sxA[k][j] * N[j]
# Counter velocity
   
return v_
def N_j(j): # Updates N based on new v
   
return N[j] + v[j] * t
for i in range(0, 25000):
    v = [v_j(
0), v_j(1), v_j(2)]
    N = [N_j(
0), N_j(1), N_j(2)]
    N_hst.append(N)
plt.plot(
range(0, len(N_hst)), [i[0] for i in N_hst])
plt.plot(
range(0, len(N_hst)), [i[1] for i in N_hst])
plt.plot(
range(0, len(N_hst)), [i[2] for i in N_hst])
plt.show()

The pattern in N was easy to identify and is application specific, but v is much more interesting as it very quickly presents itself as following Annila and Salthe’s equation of ‘v_j’ (eq.5). The code is not one hundred percent true to its origins as we have implemented a ‘counter velocity’ term. As we haven’t yet dissected our way back to the thermodynamic fundamentals, we do not yet know the relationships between respective free energies ‘A_jk‘ or conductances ‘sigma_jk‘. By implementing an opposite velocity for every determined one, we can ensure the total worth of the system is conserved at the expense of the model’s truth. Now that the number of commodities can be variable, the system can be investigated further.

# Account for increase in commodities
N = [50, 10, 40, 0, 100, 20, 50]
sig_x_A = sxA = [[
0, 1, 1, 1, 1, 1, 1],
   
[1, 0, 1, 1, 1, 1, 1],
   
[1, 1, 0, 1, 1, 1, 1],
   
[1, 1, 1, 0, 1, 1, 1],
   
[1, 1, 1, 1, 0, 1, 1],
   
[1, 1, 1, 1, 1, 0, 1],
   
[1, 1, 1, 1, 1, 1, 0]]
t =
0.0001
steps = 7500
# --------------------------------------------------------------------------
N_hst = [N]
v = [
0 for i in N]
def v_j(j):
    v_ =
0
   
for k in range(0, len(N)):
        v_ += sxA[j][k] * N[k]
        v_ += -sxA[k][j] * N[j]
   
return v_
def N_j(j):
   
return N[j] + v[j] * t
for i in range(0, steps):
    v = [v_j(i)
for i in range(0, len(v))]
    N = [N_j(i)
for i in range(0, len(N))]
    N_hst.append(N)
for l in range(0, len(N)):
    plt.plot(
range(0, len(N_hst)), [i[l] for i in N_hst])
plt.show()

No changes have been made other than the input variables. From now on, the ‘# ---…’ will mark the boundary between the two sections so identical code is not included needlessly. Again, we see that the code still works as expected when all commodities are equally valuable, they all tend to the same value.

# Investigate sxA values horizontally
N = [50, 10, 40, 0, 100, 20, 50]
sig_x_A = sxA = [[
0, 2, 2, 2, 2, 2, 2],
   
[2, 0, 2, 2, 2, 2, 2],
   
[1, 1, 0, 1, 1, 1, 1],
   
[1, 1, 1, 0, 1, 1, 1],
   
[3, 3, 3, 3, 0, 3, 3],
   
[1, 1, 1, 1, 1, 0, 1],
   
[1, 1, 1, 1, 1, 1, 0]]
t =
0.0001
steps = 7500
# --------------------------------------------------------------------------

# Investigate sxA values vertically
N = [50, 10, 40, 0, 100, 20, 50]
sig_x_A = sxA = [[
0, 1, 1, 2, 1, 1, 1],
   
[1, 0, 1, 2, 1, 1, 1],
   
[1, 1, 0, 2, 1, 1, 1],
   
[1, 1, 1, 0, 1, 1, 1],
   
[1, 1, 1, 2, 0, 1, 1],
   
[1, 1, 1, 2, 1, 0, 1],
   
[1, 1, 1, 2, 1, 1, 0]]
t =
0.0001
steps = 7500
# --------------------------------------------------------------------------

When there is a scale applied horizontally to the sxA matrix, it proportionally increases the value of the corresponding commodity. We can be confident in this, despite not having yet labelled the commodities, by the number of commodities that tend to a given value.

Using the same idea, a vertical scaling proportionally decreases the corresponding commodity’s value.

# Investigate sxA values randomly
def r():
   
return random.random()
N = [
100*r() for i in range(0, 7)]
sig_x_A = sxA = [[r()
for i in range(0, 7)],
   
[r() for i in range(0, 7)],
   
[r() for i in range(0, 7)],
   
[r() for i in range(0, 7)],
   
[r() for i in range(0, 7)],
   
[r() for i in range(0, 7)],
   
[r() for i in range(0, 7)],
]
t =
0.0001
steps = 20000
# --------------------------------------------------------------------------

Console:

sxA table 7x7

[0.492, 0.466, 0.313, 0.482, 0.753, 0.543, 0.793]
   [0.645, 0.99, 0.139, 0.713, 0.96, 0.56, 0.407]
   [0.105, 0.875, 0.126, 0.946, 0.605, 0.03, 0.86]
   [0.12, 0.791, 0.476, 0.679, 0.778, 0.694, 0.547]
   [0.518, 0.57, 0.192, 0.966, 0.086, 0.524, 0.004]
   [0.612, 0.949, 0.973, 0.759, 0.296, 0.223, 0.139]
   [0.073, 0.757, 0.52, 0.476, 0.767, 0.308, 0.864]

Introducing random values creates much more interesting systems. So far, we have only observed figures that look very alike to lnx or 1/x, but now we see commodities like the pink here, which start under-valued, become over-valued, and then settle somewhere in-between. However, Annila & Salthe did not suggest an input of a combined sxA but the ability to inform ‘sigma_jk’ and ‘A_jk’ separately. This is where, however, one could imagine an artificial intelligence could be applied, to observe an initial section of these figures, and attempt to choose a sxA ‘matrix’ which agrees and completes the system and predicts the commodities’ worth.

# Separate sig_x_A
def r():
   
return random.random()
N = [
100*r() for i in range(0, 7)]
sig = [[
0, 1, 1, 1, 1, 1, 1],
   
[1, 0, 1, 1, 1, 1, 1],
   
[1, 1, 0, 1, 1, 1, 1],
   
[1, 1, 1, 0, 1, 1, 1],
   
[1, 1, 1, 1, 0, 1, 1],
   
[1, 1, 1, 1, 1, 0, 1],
   
[1, 1, 1, 1, 1, 1, 0]]
A = [[
0, 1, 1, 1, 1, 1, 1],
   
[1, 0, 1, 1, 1, 1, 1],
   
[1, 1, 0, 1, 1, 1, 1],
   
[1, 1, 1, 0, 1, 1, 1],
   
[1, 1, 1, 1, 0, 1, 1],
   
[1, 1, 1, 1, 1, 0, 1],
   
[1, 1, 1, 1, 1, 1, 0]]
sxA = []
for j in range(0, len(sig)):
    sxA.append([])
   
for k in range(0, len(sig)):
         sxA[j].append(sig[j][k]*A[j][k])
t =
0.0001
steps = 20000
# --------------------------------------------------------------------------

We can begin to separate ‘sigma_jk’ and ‘A_jk’, but it is important we understand from the maths section how we arrive at .

# Devolve back to thermodynamic analogies
def r():
   
return random.random()
kBT =
1
n = 30 # The number of commodities
N = [r() for i in range(0, n)] # Random initial worths
sig = []
for j in range(0, n):
    sig.append([])
   
for k in range(0, n):
        sig[j].append(r())
# Random conductances
Q = []
for j in range(0, n):
    Q.append([])
   
for k in range(0, n):
        Q[j].append(r())
# Random energy pairings
t =
0.0001
steps = 7500
# --------------------------------------------------------------------------
N_hst = [N]
v = [
0 for i in N]
def v_j(j):
    v_ =
0
   
for k in range(0, len(N)):
        v_ += sxA[j][k] * N[k]
        v_ += -sxA[k][j] * N[j]
   
return v_
def N_j(j):
   
return N[j] + v[j] * t
for m in range(0, steps):
    G = [np.log(
1 / N[i]) for i in range(0, n)] # Updating the Gibbs every iteration as it depends on N
    phi = [N[i] * (np.exp(G[i] / kBT))
for i in range(0, n)] # Updating phi – energy density
    mu = [kBT * np.log(phi[i])
for i in range(0, n)] # Updating mu - potential
    A = []
   
for j in range(0, n):
        A.append([])
       
for k in range(0, n):
            A[j].append(mu[k] + Q[j][k] - mu[j])
# Free energy now has to be updated every cycle...
    sxA = []
   
for j in range(0, n):
        sxA.append([])
   
for k in range(0, n):
        sxA[j].append(sig[j][k] * A[j][k])
# ...as well as sxA=
    v = [v_j(i)
for i in range(0, n)]
    N = [N_j(i)
for i in range(0, n)]
    N_hst.append(N)
for l in range(0, n):
    plt.plot(
range(0, len(N_hst)), [i[l] for i in N_hst])
plt.show()
print(f'sxA table {len(sxA)}x{len(sxA)}')
for i in range(0, len(sxA)):
   
print([round(j * 10) / 10 for j in sxA[i]])

We now have a system that can (theoretically) handle as many commodities as we wish and balance them all in accordance with our thermoeconomic model (eq.9), (eq.10). There is likely more that could be done both physically and computationally to probe deeper into this model. The Gibbs here is determined by a relationship with log(1/N) which was a trial-and-error solution that found a stable system and produced a familiar graph. It is likely that G could be determined in another, more obviously physical, way, but it is beyond the scope of this content.

The code is easily modifiable and one is encouraged to play with it, even if just to see how the random generation affects the end result.

System requirements: Python 3 with modules: Matplotlib.pyplot , Numpy , Random.

Back to thermoeconomics

The Evolution of Physics in Finance

Completed by: Jake Wilkes, James Penston, Sam Jones and James Lundie.

Facebook icon
Instagram icon
Twitter icon

© 2022 Physics in society project

Intuit Mailchimp logo