0%

Energy Economic Codes

Posted on Edited on In ,
Symbols count in article: 21k Reading time ≈ 19 mins.

1 Markov model

1.1 Preparation and definition

  • Import modules

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    import numpy as np
    import pandas as pd
    from numpy import linalg as la

    from sklearn.metrics import mean_absolute_error, r2_score
    from sklearn.metrics import mean_squared_error

    # Plot module
    import matplotlib.pyplot as plt
    plt.rcParams['font.sans-serif'] = ['SimHei']
    # Microsoft YaHei, Times New Roman
    plt.rcParams['axes.unicode_minus'] = False

  • Plot results

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    def plot_result(test_Y, pred_Y, Title):
        ylim_min = np.round(test_Y.min()) - 1
        ylim_max = np.round(test_Y.max()) + 1
        
        # plot comparison results
        fig = plt.figure(figsize = (8, 6))
        ax0 = fig.add_subplot(111)
        ax0.scatter(test_Y, pred_Y)
        ax0.plot([ylim_min, ylim_max], [ylim_min, ylim_max], '--k')
        ax0.set_xlabel('True Target', fontsize = 15)
        ax0.set_ylabel('Target predicted', fontsize = 15)
        ax0.set_title(Title , fontsize = 16)
        ax0.text(ylim_min + 2, ylim_max - 4, r'$R^2$ = %.2f, $MAE$ = %.2f' %(r2_score(test_Y, pred_Y), 
                    mean_absolute_error(test_Y, pred_Y)), fontsize = 14)
        ax0.text(ylim_min + 2, ylim_max - 5, r'$MSE$ = %.2f, $RMSE$ = %.2f' %(mean_squared_error(test_Y, pred_Y), 
                    mean_squared_error(test_Y, pred_Y, squared=False)), fontsize = 14)
        
        ax0.set_xlim([ylim_min, ylim_max])
        ax0.set_ylim([ylim_min, ylim_max])
        plt.show()

        # Plot prediction and real values
        plt.figure(figsize = (8, 6)) # length x width
        plt.plot(pred_Y, 'r', label='prediction', lw = 0.8)
        plt.plot(test_Y, 'b', label='real', lw = 0.8)
        plt.xticks(range(19), fontsize = 15)
        plt.yticks(fontsize = 15)
        plt.title(Title, fontsize = 16)
        plt.legend(loc='best', fontsize  = 15)
        plt.show()

    绘制结果图（可选，可通过调整 plot_b 的值来选择是否绘图），包含两个:

    Fig. 1-1 Trend figure of Coal

    Fig. 1-2 Comparison figure of Coal

  • Calculate loss

    1
    2
    3
    4
    5
    6
    7
    8
    9
    def calculate_error(test_data, np_Pred_results, title):
        # Storing error results 
        R2 = r2_score(test_data, np_Pred_results)
        MAE = mean_absolute_error(test_data, np_Pred_results)
        MSE = mean_squared_error(test_data, np_Pred_results)
        RMSE = mean_squared_error(test_data, np_Pred_results, squared=False)
        
        error = [title, R2, MAE, MSE, RMSE]
        return error

    Calculate the loss: R2, MAE, MSE, RMSE

  • Definition of Markov model

    The first step: calculate and return transfer matrix (type: np.array)

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    def Markov_trans(data_stru, year_start, year_end):
        P = np.identity(4)
        K = np.array([1.0, 1., 1., 1.])
        
        year_start_data = data_stru.loc[year_start, :]
        year_end_data = data_stru.loc[year_end,:]
        
        bool_value = np.array(year_end_data < year_start_data)
        i = 0 # index
        for b in bool_value:
            if b == False: # Decrease
                K[i] =  year_end_data[i] - year_start_data[i]
            elif b == True: # Increase
                P[i, i] = year_end_data[i]/year_start_data[i]
                K[i] = 0
            i += 1
              
        KK = sum(K)
        
        for j in range(4):
            for n in range(4):
                if (j != n) & (K[j] == 0):
                    P[j,n] = (1-P[j,j])*K[n]/KK
        
        # v 为特征值    Q 为特征向量
        P = np.round(P, 3)
        v, Q = la.eig(P)
        # diag_P = np.round(np.dot(np.dot(la.inv(Q), P), Q), 3)
        V = np.diag(v)**((1)/(year_end-year_start))
        Predict_P =  np.round(np.dot(np.dot(Q, V), la.inv(Q)), 3)
        return Predict_P

    The second step: calculate and return forecasting result (type: np.array)

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    def Markov_predict(data_stru, P, year_end, year_pred):
        year_end_data = data_stru.loc[year_end,:]
        P = np.round(P, 3)
        v, Q = la.eig(P)
        # diag_P = np.round(np.dot(np.dot(la.inv(Q), P), Q), 3)
        V = np.diag(v)**(year_pred - year_end)
        Predict_P =  np.round(np.dot(np.dot(Q, V), la.inv(Q)), 3)
        np_year_end_data =np.array(year_end_data)
        Predict_energy_stru = np.around(np.dot(np_year_end_data, Predict_P), 3)
        return Predict_energy_stru

1.2 Preparation

  • Read data file

    Data.xlsx 文件需要根程序文件在同一目录

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    # Read datafile
    data_stru = pd.read_excel('./Data.xlsx', header = 0, 
                              sheet_name = 0, index_col= 0)
    num_year = len(data_stru) # the length of the dataset
    train_percent = 0.7 # Set the first 70% data as train set
    Error_results = [] # Store all error results
    plot_b = False # Choose whether plot result

    # Dataset start with 1953, end with 2017
    year_train_start, year_end = 1953, 2018
    year_train_end = year_train_start + round(len(data_stru) * 0.7) # = 1999

    test_data = np.array(data_stru.loc[year_train_end:, :]) # from 1999 to 2017

1.3 Model application

  • Avarage trans_P

    Model test based on avarage trans_P

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    trans_P = [] # store all transfer matric
    for startY in range(year_train_start, year_train_end - 1):
        for endY in range(startY, year_train_end):
            trans_temp = Markov_trans(year_start = year_train_start, 
                            year_end = year_train_end, data_stru = data_stru)
            temp = trans_temp.tolist()
            trans_P.append(temp)

    np_trans_P = np.array(trans_P) # Transfer list into numpy
    av_trans_P = np.sum(np_trans_P, 0)/np_trans_P.shape[0] # Calculate avarage

    Pred_results = [] # Using average transfer matrix to predict energy structure
    for predY in range(year_train_end, year_end):
        # Forecasting
        Pred_temp = Markov_predict(data_stru = data_stru, P = av_trans_P, 
                                year_end = year_train_end - 1, year_pred = predY)
        temp = Pred_temp.tolist() # Turn numpy.array into list
        Pred_results.append(temp) # It's convenient to append value using list
        
    np_Pred_results = np.array(Pred_results) # Turn list into np.array
    Columns = data_stru.columns
    # Plot forecasting result if needed
    if plot_b == True:
        for i in range(4):
            test_Y, pred_Y = test_data[:, i], np_Pred_results[:, i]
            Title = '{} based on average transfer matrix to predict:{}'.format(regression, 
                                        Columns[i])
            plot_result(test_Y, pred_Y, Title)
        
    # Storing error results
    title = 'average transfer matrix'  
    # Calculate 4 error
    error = calculate_error(test_data, np_Pred_results, title)
    Error_results.append(error)

  • One step and 46 step

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    for year_train_start in [1953, 1998]: # from '1953 to 1999' or '1998 to 1999'
        year_train_end, year_end = 1999, 2018
        trans_temp = Markov_trans(year_start = year_train_start, 
                                year_end = year_train_end, data_stru = data_stru)
        Pred_results = []
        for predY in range(year_train_end, year_end):
            Pred_temp = Markov_predict(data_stru = data_stru, P = trans_temp, 
                                    year_end = year_train_end - 1, year_pred = predY)
            temp = Pred_temp.tolist()
            Pred_results.append(temp)
            
        np_Pred_results = np.array(Pred_results)
        Columns = data_stru.columns
        if plot_b == True:
            for i in range(4):
                test_Y, pred_Y = test_data[:, i], np_Pred_results[:, i]
                Title = '{} based on {} and {} to predict:{}'.format(regression, 
                                        year_train_start, year_train_end, Columns[i])             
                plot_result(test_Y, pred_Y, Title)
        
        # Storing error results  
        title = '{} step'.format(year_train_end - year_train_start)   
        error = calculate_error(test_data, np_Pred_results, title)
        Error_results.append(error)

1.4 Export result to excel

1
2
3
4
np_error = np.array(Error_results)
columns = ['Transfer Matrix', 'R2', 'MAE', 'MSE', 'RMSE']
pd_error = pd.DataFrame(np_error, columns = columns, dtype = 'object')
pd_error.to_excel('Error result.xlsx', encoding = 'utf_8_sig')

1.5 All codes

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
# Markov model for energy structure prediction#
# ========================= Import modules ====================================
import numpy as np
import pandas as pd
from numpy import linalg as la

from sklearn.metrics import mean_absolute_error, r2_score
from sklearn.metrics import mean_squared_error

# Plot module
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif'] = ['SimHei']
# Microsoft YaHei, Times New Roman
plt.rcParams['axes.unicode_minus'] = False

# ======================== Preparation =======================================

# Read datafile
data_stru = pd.read_excel('./Data.xlsx', header = 0, 
                          sheet_name = 0, index_col= 0)
num_year = len(data_stru) # the length of the dataset
train_percent = 0.7 # Set the first 70% data as train set
Error_results = [] # Store all error results
plot_b = False # Choose whether plot result

# Dataset start with 1953, end with 2017
year_train_start, year_end = 1953, 2018
year_train_end = year_train_start + round(len(data_stru) * 0.7) # = 1999

test_data = np.array(data_stru.loc[year_train_end:, :]) # from 1999 to 2017

# ======================== Plot Results =====================================
def plot_result(test_Y, pred_Y, Title):
    ylim_min = np.round(test_Y.min()) - 1
    ylim_max = np.round(test_Y.max()) + 1
    
    # plot comparison results
    fig = plt.figure(figsize = (8, 6))
    ax0 = fig.add_subplot(111)
    ax0.scatter(test_Y, pred_Y)
    ax0.plot([ylim_min, ylim_max], [ylim_min, ylim_max], '--k')
    ax0.set_xlabel('True Target', fontsize = 15)
    ax0.set_ylabel('Target predicted', fontsize = 15)
    ax0.set_title(Title , fontsize = 16)
    ax0.text(ylim_min + 2, ylim_max - 4, r'$R^2$ = %.2f, $MAE$ = %.2f' %(r2_score(test_Y, pred_Y), 
                mean_absolute_error(test_Y, pred_Y)), fontsize = 14)
    ax0.text(ylim_min + 2, ylim_max - 5, r'$MSE$ = %.2f, $RMSE$ = %.2f' %(mean_squared_error(test_Y, pred_Y), 
                mean_squared_error(test_Y, pred_Y, squared=False)), fontsize = 14)
    
    ax0.set_xlim([ylim_min, ylim_max])
    ax0.set_ylim([ylim_min, ylim_max])
    plt.show()

    # Plot prediction and real values
    plt.figure(figsize = (8, 6)) # length x width
    plt.plot(pred_Y, 'r', label='prediction', lw = 0.8)
    plt.plot(test_Y, 'b', label='real', lw = 0.8)
    plt.xticks(range(19), fontsize = 15)
    plt.yticks(fontsize = 15)
    plt.title(Title, fontsize = 16)
    plt.legend(loc='best', fontsize  = 15)
    plt.show()

# Calculate the loss: R2, MAE, MSE, RMSE
def calculate_error(test_data, np_Pred_results, title):
    # Storing error results 
    R2 = r2_score(test_data, np_Pred_results)
    MAE = mean_absolute_error(test_data, np_Pred_results)
    MSE = mean_squared_error(test_data, np_Pred_results)
    RMSE = mean_squared_error(test_data, np_Pred_results, squared=False)
    
    error = [title, R2, MAE, MSE, RMSE]
    return error

# =================== Markov model definition =============================
# The first step: calculate and return transfer matrix (type: np.array)
def Markov_trans(data_stru, year_start, year_end):
    P = np.identity(4)
    K = np.array([1.0, 1., 1., 1.])
    
    year_start_data = data_stru.loc[year_start, :]
    year_end_data = data_stru.loc[year_end,:]
    
    bool_value = np.array(year_end_data < year_start_data)
    i = 0 # index
    for b in bool_value:
        if b == False: # Decrease
            K[i] =  year_end_data[i] - year_start_data[i]
        elif b == True: # Increase
            P[i, i] = year_end_data[i]/year_start_data[i]
            K[i] = 0
        i += 1
          
    KK = sum(K)
    
    for j in range(4):
        for n in range(4):
            if (j != n) & (K[j] == 0):
                P[j,n] = (1-P[j,j])*K[n]/KK
    
    # v 为特征值    Q 为特征向量
    P = np.round(P, 3)
    v, Q = la.eig(P)
    # diag_P = np.round(np.dot(np.dot(la.inv(Q), P), Q), 3)
    V = np.diag(v)**((1)/(year_end-year_start))
    Predict_P =  np.round(np.dot(np.dot(Q, V), la.inv(Q)), 3)
    return Predict_P

# The second step: calculate and return forecasting result (type: np.array)    
def Markov_predict(data_stru, P, year_end, year_pred):
    year_end_data = data_stru.loc[year_end,:]
    P = np.round(P, 3)
    v, Q = la.eig(P)
    # diag_P = np.round(np.dot(np.dot(la.inv(Q), P), Q), 3)
    V = np.diag(v)**(year_pred - year_end)
    Predict_P =  np.round(np.dot(np.dot(Q, V), la.inv(Q)), 3)
    np_year_end_data =np.array(year_end_data)
    Predict_energy_stru = np.around(np.dot(np_year_end_data, Predict_P), 3)
    return Predict_energy_stru
  
# Lable
regression = 'Markov model'

# =========== Model test based on avarage trans_P ==========================
trans_P = [] # store all transfer matric
for startY in range(year_train_start, year_train_end - 1):
    for endY in range(startY, year_train_end):
        trans_temp = Markov_trans(year_start = year_train_start, 
                        year_end = year_train_end, data_stru = data_stru)
        temp = trans_temp.tolist()
        trans_P.append(temp)

np_trans_P = np.array(trans_P) # Transfer list into numpy
av_trans_P = np.sum(np_trans_P, 0)/np_trans_P.shape[0] # Calculate avarage

Pred_results = [] # Using average transfer matrix to predict energy structure
for predY in range(year_train_end, year_end):
    # Forecasting
    Pred_temp = Markov_predict(data_stru = data_stru, P = av_trans_P, 
                            year_end = year_train_end - 1, year_pred = predY)
    temp = Pred_temp.tolist() # Turn numpy.array into list
    Pred_results.append(temp) # It's convenient to append value using list
    
np_Pred_results = np.array(Pred_results) # Turn list into np.array
Columns = data_stru.columns
# Plot forecasting result if needed
if plot_b == True:
    for i in range(4):
        test_Y, pred_Y = test_data[:, i], np_Pred_results[:, i]
        Title = '{} based on average transfer matrix to predict:{}'.format(regression, 
                                    Columns[i])
        plot_result(test_Y, pred_Y, Title)
    
# Storing error results
title = 'average transfer matrix'  
# Calculate 4 error
error = calculate_error(test_data, np_Pred_results, title)
Error_results.append(error)

# ======== Model test based on datas from 1953 to 1999 ==================
for year_train_start in [1953, 1998]: # from '1953 to 1999' or '1998 to 1999'
    year_train_end, year_end = 1999, 2018
    trans_temp = Markov_trans(year_start = year_train_start, 
                            year_end = year_train_end, data_stru = data_stru)
    Pred_results = []
    for predY in range(year_train_end, year_end):
        Pred_temp = Markov_predict(data_stru = data_stru, P = trans_temp, 
                                year_end = year_train_end - 1, year_pred = predY)
        temp = Pred_temp.tolist()
        Pred_results.append(temp)
        
    np_Pred_results = np.array(Pred_results)
    Columns = data_stru.columns
    if plot_b == True:
        for i in range(4):
            test_Y, pred_Y = test_data[:, i], np_Pred_results[:, i]
            Title = '{} based on {} and {} to predict:{}'.format(regression, 
                                    year_train_start, year_train_end, Columns[i])             
            plot_result(test_Y, pred_Y, Title)
    
    # Storing error results  
    title = '{} step'.format(year_train_end - year_train_start)   
    error = calculate_error(test_data, np_Pred_results, title)
    Error_results.append(error)

np_error = np.array(Error_results)
columns = ['Transfer Matrix', 'R2', 'MAE', 'MSE', 'RMSE']
pd_error = pd.DataFrame(np_error, columns = columns, dtype = 'object')
pd_error.to_excel('Error result.xlsx', encoding = 'utf_8_sig')

2 Traditional models

  • Ridge
  • Lasso
  • PLS
  • MLP
  • Linear SVM
  • DecisionTree

不同之处在于104 行开始的模型定义不同

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
# ========================= Import modules ====================================
import numpy as np
import pandas as pd

from sklearn.model_selection import train_test_split
from sklearn.linear_model import Ridge, Lasso
from sklearn.neural_network import MLPRegressor
from sklearn.cross_decomposition import PLSRegression
from sklearn.svm import LinearSVR
from sklearn.tree import DecisionTreeRegressor

from sklearn.metrics import mean_absolute_error, r2_score
from sklearn.metrics import mean_squared_error

# Plot module
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif'] = ['SimHei']
# Microsoft YaHei, Times New Roman
plt.rcParams['axes.unicode_minus'] = False

# ======================== Preparation =======================================

# Read datafile
data_stru = pd.read_excel('./Data.xlsx', header = 0, 
                          sheet_name = 0, index_col= 0)
data_stru = data_stru.dropna()
dataset = data_stru.values
dataset = dataset.astype('float32')

Columns = data_stru.columns # columns name
num_year = len(data_stru) # the length of the dataset
train_percent = 0.70 # Set the first 70% data as train set
test_percent = 1 - train_percent
Error_results = [] # Store all error results
input_size = 10 # Feed the previous 10 years as the input
num_features = 4 # 4 energy categories
plot_b = False # Choose whether plot result

# # Put NN calculating module into GPU if available
# device = torch.device("cuda" if torch.cuda.is_available() else "cpu")

# Feed the previous 10 years as the input and the 11th year as prediction output
def create_dataset(dataset, look_back = input_size):
    dataX, dataY = [], []
    for i in range(len(dataset) - look_back):
        a = dataset[i:(i + look_back)]
        dataX.append(a)
        dataY.append(dataset[i + look_back])
    return np.array(dataX), np.array(dataY)

# Input datasets and output datasets
data_X, data_Y = create_dataset(dataset)
ylim_min = np.round(data_Y.min()) - 1
ylim_max = np.round(data_Y.max()) + 1

# Divide training set and test set
train_X, test_X, train_Y, test_Y = train_test_split(data_X, data_Y, 
        shuffle = False, test_size = test_percent) 
# shuffle: don't shuffle data


# ======================== Plot Results =====================================
def plot_result(test_Y, pred_Y, Title):
    ylim_min = np.round(test_Y.min()) - 1
    ylim_max = np.round(test_Y.max()) + 1
    
    # plot comparison results
    fig = plt.figure(figsize = (8, 6))
    ax0 = fig.add_subplot(111)
    ax0.scatter(test_Y, pred_Y)
    ax0.plot([ylim_min, ylim_max], [ylim_min, ylim_max], '--k')
    ax0.set_xlabel('True Target', fontsize = 15)
    ax0.set_ylabel('Target predicted', fontsize = 15)
    ax0.set_title(Title , fontsize = 16)
    ax0.text(ylim_min + 2, ylim_max - 4, r'$R^2$ = %.2f, $MAE$ = %.2f' %(r2_score(test_Y, pred_Y), 
                mean_absolute_error(test_Y, pred_Y)), fontsize = 14)
    ax0.text(ylim_min + 2, ylim_max - 5, r'$MSE$ = %.2f, $RMSE$ = %.2f' %(mean_squared_error(test_Y, pred_Y), 
                mean_squared_error(test_Y, pred_Y, squared=False)), fontsize = 14)
    
    ax0.set_xlim([ylim_min, ylim_max])
    ax0.set_ylim([ylim_min, ylim_max])
    plt.show()

    # Plot prediction and real values
    plt.figure(figsize = (8, 6)) # length x width
    plt.plot(pred_Y, 'r', label='prediction', lw = 0.8)
    plt.plot(test_Y, 'b', label='real', lw = 0.8)
    plt.xticks(range(19), fontsize = 15)
    plt.yticks(fontsize = 15)
    plt.title(Title, fontsize = 16)
    plt.legend(loc='best', fontsize  = 15)
    plt.show()

# Calculate the loss: R2, MAE, MSE, RMSE
def calculate_error(test_data, np_Pred_results, title):
    # Storing error results
    R2 = float('{:.5f}'.format(r2_score(test_data, np_Pred_results)))
    MAE = float('{:.5f}'.format(mean_absolute_error(test_data, np_Pred_results)))
    MSE = float('{:.5f}'.format(mean_squared_error(test_data, np_Pred_results)))
    RMSE =float('{:.5f}'.format( mean_squared_error(test_data, np_Pred_results, squared=False)))
    error = [title, R2, MAE, MSE, RMSE]
    return error

# ================= Ridge Regression =============================
# Ridge model
def predict_model(model_train_X, model_train_Y):
    regr = Ridge()
    regr.fit(model_train_X, model_train_Y)
    model_pred_Y = regr.predict(model_test_X)
    regression_name = 'Ridge Regression'
    return model_pred_Y, regression_name

# # Lasso model
# def predict_model(model_train_X, model_train_Y):
#     regr = Lasso()
#     regr.fit(model_train_X, model_train_Y)
#     model_pred_Y = regr.predict(model_test_X)
#     regression_name = 'Lasso Model'
#     return model_pred_Y, regression_name

# # PLS model
# def predict_model(model_train_X, model_train_Y):
#     regr = PLSRegression(n_components=1)
#     regr.fit(model_train_X, model_train_Y)
#     model_pred_Y = regr.predict(model_test_X)
#     regression_name = 'PLS Regression'
#     return model_pred_Y, regression_name

# ================= MLP Regression =============================
# # MLP model
# def predict_model(model_train_X, model_train_Y):
#     regr = MLPRegressor(random_state = 1, max_iter = 500)
#     regr.fit(model_train_X, model_train_Y)
#     model_pred_Y = regr.predict(model_test_X)
#     regression_name = 'MLP Regression'
#     return model_pred_Y, regression_name

# ================= Linear SVM regression =============================
# # SVC model
# def predict_model(model_train_X, model_train_Y):
#     regr = LinearSVR(random_state= 0, tol= 1e-5)
#     regr.fit(model_train_X, model_train_Y)
#     model_pred_Y = regr.predict(model_test_X)
#     regression_name = 'Linear SVR Regression'
#     return model_pred_Y, regression_name

# ================= DecisionTree Regressor =============================
# # DecisionTree model
# def predict_model(model_train_X, model_train_Y):
#     regr = DecisionTreeRegressor(random_state=0)
#     regr.fit(model_train_X, model_train_Y)
#     model_pred_Y = regr.predict(model_test_X)
#     regression_name = 'Decision-Tree Regression'
#     return model_pred_Y, regression_name

# Forecasting
pred_result = [] # to store prediction results
for i in range(4):
    if i < 3:
        model_train_X = train_X[:, :, i].reshape(-1, input_size)
        model_test_X = test_X[:, :, i].reshape(-1, input_size)
        model_train_Y = train_Y[:, i].reshape(-1)
        model_test_Y = test_Y[:, i].reshape(-1)
        model_pred_Y, regression = predict_model(model_train_X = model_train_X,
                    model_train_Y = model_train_Y)
        pred_result.append(model_pred_Y.tolist())
    elif i == 3:
        # Calculate the last energy consumption percent
        np_pred_result = np.array(pred_result).T
        if len(np_pred_result.shape) == 3: # Model output is 3-dimension
            model_pred_Y = 100 - np_pred_result.sum(2).reshape(-1, 1) # PLS model
            np_pred_result = np.concatenate((np_pred_result[0], model_pred_Y), axis = 1) # for PLS model
        elif len(np_pred_result.shape) == 2: # Model output is 2-dimension
            model_pred_Y = 100 - np_pred_result.sum(1).reshape(-1, 1) # Lasso model 
            np_pred_result = np.concatenate((np_pred_result, model_pred_Y), axis = 1) # for Lasso model
    Title = 'Using {} to predict:{}'.format(regression, Columns[i])
    if plot_b == True:
        plot_result(model_test_Y, model_pred_Y, Title) 
        
# Updata Error table
error = calculate_error(test_Y, np_pred_result, regression)
np_error = np.array(error).reshape(1, -1)
columns = ['Transfer Matrix', 'R2', 'MAE', 'MSE', 'RMSE']
pd_error = pd.DataFrame(np_error, columns = columns)

error_excel =  pd.read_excel('./Error result.xlsx', header = 0, sheet_name = 0, index_col= None)
res = pd.concat([error_excel, pd_error], axis = 0, ignore_index = True)

# Updata predict table
model_col = []
for c in Columns:
    temp = '{}_{}'.format(c, regression.split(' ')[0])
    model_col.append(temp)
pd_pred_result = pd.DataFrame(np_pred_result, columns = model_col, index = range(2001, 2018))
pred_excel =  pd.read_excel('./Predict results.xlsx', header = 0, sheet_name = 0, index_col= 0)
pred_res = pd.concat([pred_excel, pd_pred_result], axis = 1)

# Store updated data
res.to_excel('./Error result.xlsx', encoding = 'utf_8_sig', header = True, index = False)
pred_res.to_excel('./Predict results.xlsx', encoding = 'utf_8_sig',  header = True, index = True)

3 Linear NN model

-------------This blog is over! Thanks for your reading-------------