Python implements threshold regression

The basic idea of the Threshold Regressive Model (TR model or TRM) is to use the control function of the threshold variable. When the predictor data is given, the threshold value of the threshold variable first determines the control function according to the threshold value of the threshold variable. Use different forecasting equations to try to explain various phenomena like jumps and sudden changes. Its essence is to classify the forecasting problems according to the value of the state space, and use the piecewise linear regression model to describe the overall nonlinear forecasting problem.

The modeling step of multiple threshold regression is the process of confirming the threshold variables, calibrating the threshold number L, the threshold value and the regression coefficient. For the convenience of calculation, two divisions (ie L=2) are used to illustrate the modeling steps of the model.

The basic steps are as follows (with code):

1. Read the data and calculate the correlation coefficient matrix between the forecast object and the forecast factor.

Data read
# Use pandas to read csv, the data read is a DataFrame object
data = pd.read_csv('jl.csv')
# Convert DataFrame object to array,The first column of the array is the data number, the last column is the forecast object, and the middle columns are the forecast factors
data= data.values.copy()
# print(data)
# Calculate the cross-correlation coefficient, the parameters are the predictor sequence and the lag time k
def get_regre_coef(X,Y,k):
 S_xy=0
 S_xx=0
 S_yy=0
 # Calculate the mean value of forecast factors and forecast objects
 X_mean = np.mean(X)
 Y_mean = np.mean(Y)for i inrange(len(X)-k):
 S_xy +=(X[i]- X_mean)*(Y[i+k]- Y_mean)for i inrange(len(X)):
 S_xx +=pow(X[i]- X_mean,2)
 S_yy +=pow(Y[i]- Y_mean,2)return S_xy/pow(S_xx*S_yy,0.5)
# Calculate the correlation coefficient matrix
def regre_coef_matrix(data):
 row=data.shape[1]#Number of columns
 r_matrix=np.ones((1,row-2))
 # print(row)for i inrange(1,row-1):
 r_matrix[0,i-1]=get_regre_coef(data[:,i],data[:,row-1],1)#Lag time is 1 return r_matrix
r_matrix=regre_coef_matrix(data)
# print(r_matrix)
### Output###
#[[0.0489790.078299890.190057050.275012090.28604638]]

2. The correlation coefficients are sorted, and the factor with the largest correlation coefficient is used as the threshold element.

# Sort the correlation coefficient and find the largest correlation coefficient as the threshold element
def get_menxiannum(r_matrix):
 row=r_matrix.shape[1]#Number of columns
 for i inrange(row):if r_matrix.max()==r_matrix[0,i]:return i+1return-1
m=get_menxiannum(r_matrix)
# print(m)
## Output##The fifth factor has the largest correlation coefficient
#5

3. Reorder the data according to the selected threshold element factor.

# Sort the factor sequence according to the threshold element,m is the serial number of the threshold variable
def resort_bymenxian(data,m):
 data=data.tolist()#Convert to list
 data.sort(key=lambda x: x[m])#List according to m+1 column to sort(Ascending)
 data=np.array(data)return data
data=resort_bymenxian(data,m)#Get the sorted sequence array

4. The sorted sequence is divided into two segments according to the threshold element, the first segment is 1 data, the second segment is n-1 (n is the sample capacity) data; the second segment is 2 data, The second segment of n-2 data, one-time analogy, respectively calculate the divided F statistic and select the segmentation point of the threshold element corresponding to the largest statistic as the threshold.

def get_var(x):return x.std()**2* x.size #Calculate the total variance
# Calculation of statistic F,The input data is the forecast object data sorted according to the threshold element
def get_F(Y):
 col=Y.shape[0]#Number of rows, sample size
 FF=np.ones((1,col-1))#Store statistics of different split points
 V=get_var(Y)#Calculate the total variance
 for i inrange(1,col):#1 to col-1
 S=get_var(Y[0:i])+get_var(Y[i:col])#Calculate the sum of variance within two segments
 F=(V-S)*(col-2)/S
 FF[0,i-1]=F#At this step, it is necessary to judge whether the F test is passed, and the F statistic is retained only after passing
 return FF
y=data[:,data.shape[1]-1]
FF=get_F(y)
def get_index(FF,element):#Get the index of the first occurrence of element in a one-dimensional array FF
 i=-1for item in FF.flat:
 i+=1if item==element:return i
f_index=get_index(FF,np.max(FF))#Get the largest index of the statistic F
# print(data[f_index,m-1])#The threshold element is the fifth factor, and the threshold value is 121 when substituted into the index

5. The data sequence is divided into two segments with the threshold as the dividing point, and multiple linear regression is performed respectively. Here, the linear regression module in the sklearn.linear_model module is used. Substituting the predictive factors to calculate the predicted values of the two segments respectively.

# Divide the new data sequence into two parts with the threshold as the dividing point, and perform multiple regression calculations respectively
def data_excision(data,f_index):
 f_index=f_index+1
 data1=data[0:f_index,:]
 data2=data[f_index:data.shape[0],:]return data1,data2
data1,data2=data_excision(data,f_index)
# First paragraph
def get_XY(data):
 # Array slice to assign values to variables
 Y = data[:, data.shape[1]-1] #The forecast object is in the last column
 X = data[:,1:data.shape[1]-1]#Predictor from the second column to the penultimate column
 return X, Y
X,Y=get_XY(data1)
regs=LinearRegression()
regs.fit(X,Y)
# print('First paragraph')
# print(regs.coef_)#Output regression coefficient
# print(regs.score(X,Y))#Output correlation coefficient
# Calculate predicted value
Y1=regs.predict(X)
# print('Second paragraph')
X,Y=get_XY(data2)
regs.fit(X,Y)
# print(regs.coef_)#Output regression coefficient
# print(regs.score(X,Y))#Output correlation coefficient
# Calculate predicted value
Y2=regs.predict(X)
Y=np.column_stack((data[:,0],np.hstack((Y1,Y2)))).copy()
Y=np.column_stack((Y,data[:,data.shape[1]-1]))
Y=resort_bymenxian(Y,0)

6. The predicted value and the actual value are re-sorted according to the year number, the order is restored, and the matplotlib module is used to make a comparison chart between the predicted value and the actual value.

# Recovery order
Y=resort_bymenxian(Y,0)
# print(Y.shape)
# Visualization of forecast results
plt.plot(Y[:,0],Y[:,1],'b--',Y[:,0],Y[:,2],'g')
plt.title('Comparison of predicted and measured values',fontsize=20,fontname='Times New Roman')#Add title
plt.xlabel('Years',color='gray')#Add x-axis label
plt.ylabel('Average traffic in December',color='gray')#Add y-axis label
plt.legend(['Predicted values','Measured values'])#Add legend
plt.show()

Result graph:

Data used: Quoted from "Modern Medium and Long-term Hydrological Forecasting Methods and Applications" by Zhang Shiming, Tang Chengyou's official article

The above python threshold regression method is all the content shared by the editor, I hope to give you a reference.