Leontief Model Simulation with Python Part 2: Experimental Attempt of Eigenvalue Problem when a matrix is singular i.e. non-invertible

 

This piece of my experimental attempt is based on the paper Singularity in the Discrete Dynamic Leontief Model by István Ábel1, Imre Dobos https://pp.bme.hu/so/article/view/8432/7719 .  This paper was accidentally found while searching for "Matrix Pencil for financial analysis".  Yet, this method is different from the Matrix Pencil (MP) method used in my current researching topic.  Nonetheless, this mathematical modelling applying the eigenvalue problem for the cross sectional econometric model is interesting enough to entice me. 

Honestly speaking, I have not fully understood their method behind because it is really off-topic from my current time series modelling. At the same time, I would like to keep this topic for my near future reference as I will be able to use my Linear Algebra application gained from my current research project. At least, I have attempted solving their mathematical modelling while learning by imitation. Therefore, it is likely to contain some error which I have to spend a sufficient amount of time for investigating and learning this realm of research. 

In terms of this example, the capital coefficient matrix C_i,i is singular containing a row containing zeros i.e. non-inversible. Therefore, the method of the eigenvalue problem is applied instead of using C_i,i^(-1).  Lambda (Greek letter) denotes the eigenvalue of the following equation. The error margin is kept below 1% or 5% at most.

My Python codes are displayed below:

 

# importing necessary tools

import matplotlib.pyplot as plt

import networkx as nx

import numpy as np

import pandas as pd

import random

import math

import statistics

import networkx as nx

import scipy.linalg

from scipy.stats import qmc

RefList=[

' https://en.wikipedia.org/wiki/Input%E2%80%93output_model'

,

' https://www.youtube.com/watch?v=KmVfmISjayA&t=134s'

,

' https://www.youtube.com/watch?v=z_9HwKet8G0&t=301s'

,

' https://www.sciencedirect.com/science/article/pii/S0895717710001093'

,

' https://pp.bme.hu/so/article/view/8432 https://pp.bme.hu/so/article/view/8432/7719'

]

print('Referring to the following articles and YouTube videos: \n')

for k in range(len(RefList)):

    print(RefList[k])

print(' \n')

print('x_i,t: the vector of output levels \n')

print('d_i,t: the vector of final demands (excluding investment) \n')

print('L: the Leontief input–output matrix \n')

print('C: the capital coefficient matrix \n')

print('C: The following equation is used when these matrices are non-singular i.e. inversible. \n')

# The following follows the example shown in the paper

# Index: Industrial Sectors

Sctrs=['Sec1','Sec2','Sec3']

print(f'There are {len(Sctrs)} industrial sectors. \n')

L = np.array([[0.3,0.3,0.3],[0.4,0.1,0.5],[0.3,0.5,0.2]])

Pdf_L = pd.DataFrame(data = L,index = Sctrs,columns = Sctrs)

display('The input–output matrix: L =',Pdf_L)

C = np.array([[0.3,0.4,0.45],[0,0,0],[0.6,0.8,0.9]])

Pdf_C = pd.DataFrame(data = C,index = Sctrs,columns = Sctrs)

display('The capital coefficient matrix: C =',Pdf_C)

# Identity Matrix

I=np.identity(len(C))

Pdf_I = pd.DataFrame(data = I,index = Sctrs,columns = Sctrs)

display(f'I: {len(C)} by {len(C[0])} identity matrix',Pdf_I)

print(' \n')

print('Leontief formula: x_i,t = L x_i,t + C _i,i (x_i,t+1 - x_i,t) \n')

print('Then, it is converted to:  x_i,t+1 = C_i,i^-1 {(I_i - L_i,i + C_i,i) x_i,t}    \n')

print('On the other hand, the capital coefficient matrix C_i,i is singular containing a row containing zeros i.e. non-inversible. \n ')

print('Therefore, the method of the eigenvalue problem is applied instead of using C_i,i^(-1) as follows. \n ')

print(f"Lambda (Greek letter) denotes the eigenvalue of the following equation (Ref. {RefList[4]} : ")

print('1. x_i,t+1 C_i,i = (I_i - L_i,i + C_i,i) x_i,t } \n')

print('2. Lambda x_i,t C_i,i = (I_i - L_i,i + C_i,i) x_i,t } \n')

print('3. Lambda C_i,i = (I_i - L_i,i + C_i,i) } \n')

# Eigenvalue finding

E= scipy.linalg.eig(C, (I - L + C))

for k in range(len(E[0])):

    print(round(E[0][k].real,18))

# Using the positive real number eigenvalue

EigVal=E[0].real

Where=(np.where(EigVal==max(EigVal)))[0][0]

Lambda=EigVal[Where]

print(f"The maximum eigenvalue used is Lambda= {Lambda} .")

# Corresponding Engenvector denoted as Rho (Greek letter)

ind=np.where(E[0] == Lambda)

ind[0][0]

Rho=E[1][:,ind[0][0]]

print(' \n')

# Let's denote the initial values of production of 3 sectors

x_0=np.array([392999.32,422999.27,446999.23])

Pdf_x_0 = pd.DataFrame(data = x_0,index = Sctrs)

display("Let's denote the initial values of production of 3 sectors: x_i,0 =",Pdf_x_0)

x_1=x_0/Lambda

Pdf_x_1 = pd.DataFrame(data = x_1,index = Sctrs)

display("Then, x_i,1 = x_i,0 / Lambda  =",Pdf_x_1)

Right=np.dot(L,x_0)+np.dot(C,x_1-x_0)

print(f"Then, the right hand side (I_i - L_i,i + C_i,i) x_i,t is {Right} , \n")

Left=x_0

print(f"and the left hand side (I_i - L_i,i + C_i,i) x_i,t is {Left}. \n")

print(f"The percentage difference between the left-hand side and the right-hand side is {(Left-Right)/Right}. \n")