|
| 1 | +"""Implementation of Lindner (2012) in Python with NumPy and Pandas. |
| 2 | +
|
| 3 | +Lindner, Sören, Julien Legault, and Dabo Guan. 2012. |
| 4 | +‘Disaggregating Input–Output Models with Incomplete Information’. |
| 5 | +Economic Systems Research 24 (4): 329–47. |
| 6 | +https://doi.org/10.1080/09535314.2012.689954. |
| 7 | +
|
| 8 | +The comments in this script contain the Matlab code given in the supplementary |
| 9 | +material 'cesr_a_689954_sup_27358897.docx' of Lindner (2012). |
| 10 | +
|
| 11 | +Source (accessed 06.12.2022): |
| 12 | +https://www.tandfonline.com/doi/suppl/10.1080/09535314.2012.689954 |
| 13 | +
|
| 14 | +The script contains the generation of random numbers. A random vector is |
| 15 | +generated in line 90 of the Matlab script: |
| 16 | +
|
| 17 | + `base(p,:) = rand(1,Nv)` |
| 18 | +
|
| 19 | +For verification purposes, `np.random.seed(1337)` (Python) and |
| 20 | +`rand('twister', 1337)` (Matlab) was applied. |
| 21 | +
|
| 22 | +""" |
| 23 | + |
| 24 | +import numpy as np |
| 25 | +import pandas as pd |
| 26 | + |
| 27 | +from tqdm import tqdm |
| 28 | + |
| 29 | +if False: # !!! |
| 30 | + |
| 31 | + # Switch flag for verification |
| 32 | + # Matlab equivalent: `rand('twister', 1337)` |
| 33 | + # Source: https://stackoverflow.com/a/20202330/5696601 |
| 34 | + |
| 35 | + np.random.seed(1337) |
| 36 | + |
| 37 | +# %% Loading IO data |
| 38 | + |
| 39 | +flows = pd.read_csv( |
| 40 | + # Input–output table of China (2007), in billion RMB |
| 41 | + './input/io-table-cn-2007-flows.csv', |
| 42 | + header=None |
| 43 | + ) |
| 44 | + |
| 45 | +flows_idx = pd.read_csv( |
| 46 | + './input/io-table-cn-2007-flows-idx.csv' |
| 47 | + ) |
| 48 | + |
| 49 | +flows.columns = pd.MultiIndex.from_frame(flows_idx) |
| 50 | +flows.index = pd.MultiIndex.from_frame(flows_idx.iloc[:12, :]) |
| 51 | + |
| 52 | +# Vector of final demand |
| 53 | +f = flows.loc[:, ('Final demand', 'FD')] |
| 54 | + |
| 55 | +# Vector of intermediate demand |
| 56 | +id = flows.loc[:, ('Intermediate demand', 'ID')] |
| 57 | + |
| 58 | +# Vector of total outputs |
| 59 | +x = f + id |
| 60 | + |
| 61 | +# Exchange matrix |
| 62 | +Z = flows.loc[ |
| 63 | + # Rows |
| 64 | + :, |
| 65 | + # Cols |
| 66 | + (~flows.columns.get_level_values('Cat') |
| 67 | + .isin(['ID', 'FD', 'TO'])) |
| 68 | + ] |
| 69 | + |
| 70 | +del flows_idx |
| 71 | + |
| 72 | +temp = Z.shape # Size of IO table |
| 73 | + |
| 74 | +N = temp[0] - 1 # Number of common sectors |
| 75 | + |
| 76 | +A = np.divide(Z, x) # ggregated technical coefficient matrix |
| 77 | + |
| 78 | +x_common = x[:-1] # Vector of total outputs for common sectors |
| 79 | + |
| 80 | +f_common = f[:-1] # Vector of final demand for common sectors |
| 81 | + |
| 82 | +# Note: The last sector of the table is disaggregated, |
| 83 | +# i.e. the electricity sector |
| 84 | + |
| 85 | +x_elec = x[-1] # Total output of the disaggregated sector |
| 86 | + |
| 87 | +f_elec = f[-1] # Final demand of the disaggregated sector |
| 88 | + |
| 89 | +# %% Newly formed sectors from the electricity sector |
| 90 | + |
| 91 | +# New sector weights |
| 92 | +w = pd.read_csv( |
| 93 | + './input/io-table-cn-2007-w.csv', |
| 94 | + header=None |
| 95 | + ) |
| 96 | + |
| 97 | +w = w.values.flatten() |
| 98 | + |
| 99 | +w_idx = pd.read_csv( |
| 100 | + './input/io-table-cn-2007-w-idx.csv' |
| 101 | + ) |
| 102 | + |
| 103 | +n = len(w) # Number of new sectors |
| 104 | + |
| 105 | +# Total number of sectors for the disaggregated IO table |
| 106 | +N_tot = N + n |
| 107 | + |
| 108 | +# Vector of new total sector outputs |
| 109 | +x_new = w*x_elec/1000 |
| 110 | + |
| 111 | +# Vector of disaggregated economy sector total outputs |
| 112 | +xs = np.concatenate((x_common, x_new)) |
| 113 | + |
| 114 | +f_new = w*f_elec # # Final demand of new sectors |
| 115 | + |
| 116 | +# %% Building the constraint matrix C |
| 117 | + |
| 118 | +Nv = n * N_tot + n # Number of variables |
| 119 | + |
| 120 | +Nc = N + n + 1 # Number of constraints |
| 121 | + |
| 122 | +# Vector of constraint constants |
| 123 | +q = pd.concat( |
| 124 | + [A.iloc[N, :], |
| 125 | + pd.Series(w, index=pd.MultiIndex.from_frame(w_idx))] |
| 126 | + ) |
| 127 | + |
| 128 | +# Matrix of constraints |
| 129 | +C = np.zeros((Nc, Nv)) |
| 130 | + |
| 131 | +# %%% Common sectors constraints |
| 132 | + |
| 133 | +C11 = np.zeros((N, N*n)) |
| 134 | + |
| 135 | +for ii in range(N): |
| 136 | + col_indices = range(n*(ii), n*ii+n) |
| 137 | + C11[ii, col_indices] = np.ones((1, n)) |
| 138 | + |
| 139 | +C[:N, :N*n] = C11 |
| 140 | + |
| 141 | +# %%% New sectors constraints |
| 142 | + |
| 143 | +C22 = np.zeros((1, n**2)) |
| 144 | + |
| 145 | +for ii in range(0, n): |
| 146 | + col_indices = range(n*(ii), n*ii+n) |
| 147 | + C22[0, col_indices] = w[ii]*np.ones((1, n)) |
| 148 | + |
| 149 | +C[N, N*n:N*n+n**2] = C22 |
| 150 | + |
| 151 | +# %%% Final demand constraints |
| 152 | + |
| 153 | +C31 = np.zeros((n, N*n)) |
| 154 | + |
| 155 | +for ii in range(N): |
| 156 | + col_indices = range(n*(ii-1)+3, n*ii+3) |
| 157 | + C31[:n, col_indices] = (x_common[ii]/x_elec)*np.eye(n) |
| 158 | + |
| 159 | +C32 = np.zeros((n, n**2)) |
| 160 | + |
| 161 | +for ii in range(0, n): |
| 162 | + col_indices = range(n*(ii-1)+3, n*ii+3) |
| 163 | + C32[:n, col_indices] = w[ii]*np.eye(n) |
| 164 | + |
| 165 | +C[N+1:, :N*n] = C31 |
| 166 | +C[N+1:, N*n:N*n+n**2] = C32 |
| 167 | +C[N+1:, N*n+n**2:] = np.eye(n) |
| 168 | + |
| 169 | +# %% Building the initial estimate y0 |
| 170 | + |
| 171 | +# Technical coefficient matrix of the initial estimate |
| 172 | +As_y0 = np.zeros((N_tot, N_tot)) |
| 173 | + |
| 174 | +# Common/Common part |
| 175 | +As_y0[:N, :N] = A.iloc[:N, :N] |
| 176 | + |
| 177 | +# Common/New part |
| 178 | +As_y0[:N, N:N_tot] = np.repeat(A.iloc[:N, N].to_numpy(), n).reshape(N, n) |
| 179 | + |
| 180 | +# New/Common part |
| 181 | +As_y0[N:N_tot, :N] = ( |
| 182 | + np.multiply(w, A.iloc[N, :N].to_numpy().repeat(n).reshape(N, n)).T |
| 183 | + ) |
| 184 | + |
| 185 | +# New/New part |
| 186 | +As_y0[N:N_tot, N:N_tot] = np.multiply( |
| 187 | + A.iloc[N, N], |
| 188 | + np.repeat(w, n).reshape(n, n) |
| 189 | + ) |
| 190 | + |
| 191 | +# %% Generating the orthogonal distinguishing matrix |
| 192 | + |
| 193 | +# %%% Making the constraint matrix orthogonal |
| 194 | + |
| 195 | +C_orth = C.copy() |
| 196 | + |
| 197 | +for c in tqdm(range(Nc), desc='Orthogonalize constraint matrix'): |
| 198 | + for i in range(c): |
| 199 | + |
| 200 | + # Orthogonal projection |
| 201 | + C_orth[c, :] = ( |
| 202 | + C_orth[c, :] |
| 203 | + - np.dot(C_orth[c, :], C_orth[i, :]) |
| 204 | + / np.linalg.norm(C_orth[i, :])**2 * C_orth[i, :] |
| 205 | + ) |
| 206 | + |
| 207 | +# %%% Gram-Schmidt algorithm |
| 208 | + |
| 209 | +base = np.zeros((Nv, Nv)) # Orthogonal base containing C_orth and D |
| 210 | +base[:Nc, :] = C_orth.copy() |
| 211 | + |
| 212 | +for p in tqdm(range(Nc, Nv), desc='Gram-Schmidt algorithm'): |
| 213 | + |
| 214 | + # Generate random vector |
| 215 | + base[p, :] = np.random.rand(1, Nv) |
| 216 | + |
| 217 | + for i in range(p): |
| 218 | + |
| 219 | + # Orthogonal projection on previous vectors |
| 220 | + base[p, :] -= ( |
| 221 | + np.dot(base[p, :], base[i, :]) |
| 222 | + / np.linalg.norm(base[i, :])**2 |
| 223 | + * base[i, :] |
| 224 | + ) |
| 225 | + |
| 226 | + # Normalizing |
| 227 | + base[p, :] /= np.linalg.norm(base[p, :]) |
| 228 | + |
| 229 | +# Retrieving the distinguishing matrix from the orthogonal base |
| 230 | +D = base[Nc:, :].T |
0 commit comments