-
Notifications
You must be signed in to change notification settings - Fork 0
/
cal_starma.py
executable file
·140 lines (119 loc) · 4.92 KB
/
cal_starma.py
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
from pySTARMA import stacf_stpacf,starma_model,utils
import pandas as pd
import numpy as np
import matplotlib.pylab as plt
from matplotlib.pylab import rcParams
import csv
#rcParams设定好画布的大小
rcParams['figure.figsize']=15,6
file='f://biyedata//Guangdong_AMJ.csv'
tsd_f='f://biyedata//ts_diff1.csv'
st_acf='f://biyedata//st_acf.csv'
st_pacf='f://biyedata//stpacf.csv'
st_predict='f://biyedata//stpredict2.csv'
dateparse = lambda dates: pd.datetime.strptime(dates, '%Y/%m/%d')
#---其中parse_dates 表明选择数据中的哪个column作为date-time信息，
#---index_col 告诉pandas以哪个column作为 index
#--- date_parser 使用一个function(本文用lambda表达式代替)，使一个string转换为一个datetime变量
#时空数据列
data=pd.read_csv(file,encoding='gbk',parse_dates=['riqi'],index_col=['riqi'],date_parser=dateparse)
#空间邻接矩阵
'''
weight=[[0,0,0.3333,0,0.3333,0,0,0,0,0.3333,0,0,0,0],
[0,0,0.1667,0,0,0.1667,0.1667,0,0,0,0,0.1667,0.1667,0.1667],
[0.1667,0.1667,0,0,0.1667,0,0.1667,0,0.1667,0,0,0,0,0.1667],
[0,0,0,0,0.3333,0,0,0.3333,0.3333,0,0,0,0,0],
[0.1667,0,0.1667,0.1667,0,0,0,0.1667,0.1667,0.1667,0,0,0,0],
[0,0.25,0,0,0,0,0,0,0,0,0.25,0.25,0,0.25],
[0,0.25,0.25,0,0,0,0,0,0.25,0,0,0,0.25,0],
[0,0,0,0.5,0.5,0,0,0,0,0,0,0,0,0],
[0,0,0.25,0.25,0.25,0,0.25,0,0,0,0,0,0,0],
[0.5,0,0,0,0.5,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0.5,0,0,0,0,0,0.5,0,0],
[0,0.25,0,0,0,0.25,0,0,0,0,0.25,0,0.25,0],
[0,0.3333,0,0,0,0,0.3333,0,0,0,0,0.3333,0,0],
[0,0.3333,0.3333,0,0,0.3333,0,0,0,0,0,0,0,0]]
'''
weight=[
[1,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,1,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,1,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,1,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,1,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,1,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,1,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,1,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,1,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,1,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,1,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,1,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,1,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,1]
]
weight1=[
[0.1111111,0.1111111,0.1111111,0.1111111,0.1111111,0.1111111,0,0,0.1111111,0,0,0.1111111,0.1111111,0.1111111],
[0.1111111,0,0.1111111,0,0.1111111,0.1111111,0.1111111,0,0,0.1111111,0,0.1111111,0.1111111,0.1111111],
[0.1111111,0.1111111,0,0,0.1111111,0.1111111,0.1111111,0,0,0.1111111,0,0.1111111,0.1111111,0.1111111],
[0.25,0,0,0,0.25,0,0,0.25,0,0.25,0,0,0,0],
[0.125,0.125,0.125,0.125,0,0.125,0,0.125,0,0.125,0,0,0,0.125],
[0.1,0.1,0.1,0,0.1,0,0.1,0,0,0.1,0.1,0.1,0.1,0.1],
[0.125,0.125,0.125,0,0,0.125,0,0,0.125,0,0,0.125,0.125,0.125],
[0,0,0,0.3333333,0.3333333,0,0,0,0,0.3333333,0,0,0,0],
[0,0,0,0,0,0,1,0,0,0,0,0,0,0],
[0.125,0.125,0.125,0.125,0.125,0.125,0,0.125,0,0,0,0,0,0.125],
[0,0,0,0,0,0.3333333,0,0,0,0,0,0.3333333,0.3333333,0],
[0,0.1428571,0.1428571,0,0,0.1428571,0.1428571,0,0,0,0.1428571,0,0.1428571,0.1428571],
[0.125,0.125,0.125,0,0,0.125,0.125,0,0,0,0.125,0.125,0,0.125],
[0.1111111,0.1111111,0.1111111,0,0.1111111,0.1111111,0.1111111,0,0,0.1111111,0,0.1111111,0.1111111,0]]
weight2=[
[0,0,0,0,0,0,0.25,0.25,0,0.25,0.25,0,0,0],
[0,0,0,0.25,0,0,0,0.25,0.25,0,0.25,0,0,0],
[0,0,0,0.25,0,0,0,0.25,0.25,0,0.25,0,0,0],
[0,0.1666667,0.1666667,0,0,0.1666667,0.1666667,0,0,0,0,0,0.1666667,0.1666667],
[0,0,0,0,0,0,0.25,0,0,0,0.25,0.25,0.25,0],
[0,0,0,0.3333333,0,0,0,0.3333333,0.3333333,0,0,0,0,0],
[0,0,0,0.25,0.25,0,0,0,0,0.25,0.25,0,0,0],
[0.2,0.2,0.2,0,0,0.2,0,0,0,0,0,0,0,0.2],
[0.1428571,0.1428571,0.1428571,0,0,0.1428571,0,0,0,0,0,0.1428571,0.1428571,0.1428571],
[0,0,0,0,0,0,0.25,0,0,0,0.25,0.25,0.25,0],
[0.1428571,0.1428571,0.1428571,0,0.1428571,0,0.1428571,0,0,0.1428571,0,0,0,0.1428571],
[0.25,0,0,0,0.25,0,0,0,0.25,0.25,0,0,0,0],
[0,0,0,0.25,0.25,0,0,0,0.25,0.25,0,0,0,0],
[0,0,0,0.25,0,0,0,0.25,0.25,0,0.25,0,0,0]]
w_slag=[weight,weight1]
a=w_slag[0]
ts=np.array(data)
w=np.array(w_slag)
#一次差分
ts_diff=utils.set_stationary(ts,[1])
#自相关函数
acf=stacf_stpacf.Stacf(ts_diff,w,5)
acf_val=acf.estimate()
#偏自相关函数
pacf=stacf_stpacf.Stpacf(ts_diff,w,5)
pacf_val=pacf.estimate()
#STARMA
#starma=starma_model.STARMA(3,1,ts_diff,w)
starma=starma_model.STARMA(1,0,ts_diff,w)
starma.fit()
starma.print_results()
print(starma.get_model())
#print(starma.predict(ts,3))
result=starma.predict(ts_diff,5)
#筛选出新的csv文件
with open(tsd_f,'w',newline='') as tsd_f:
writer=csv.writer(tsd_f)
for row in ts_diff:
writer.writerow(row)
with open(st_acf,'w',newline='')as acf_f:
writer=csv.writer(acf_f)
for row in acf_val:
writer.writerow(row)
with open(st_pacf,'w',newline='')as pacf_f:
writer=csv.writer((pacf_f))
for row in pacf_val:
writer.writerow(row)
with open(st_predict,'w',newline='')as predict_f:
writer=csv.writer((predict_f))
for row in result:
writer.writerow(row)