用數學計算日馳何時崩盤!

目前我持有日馳的股票,前段時間噴高高,但最近不停漲停!跌停!搞的我心力交瘁。於是今天決定用一個數學模型來預測可能的崩盤時機。
LPPL是我之前在ricequant看到的泡沫預測模型,對於那些極右側交易者(追漲)滿有幫助的,然而這個模型只考慮股價,不考慮財報、月報等,所以比較適合用在大盤,個股可能比較沒辦法。但我還是很想知道日馳是否要崩盤了!?XD

預防針

這邊還是要提醒大家,這個只是個有趣的小實驗,但我真的不確定這個很可靠,請不要依賴LPPL買賣股票,可以當作一個參考,追蹤一下日馳的股價是否跟LPPL模型顯示的一樣。我個人是持非常懷疑的態度!請大家不要因為這篇文章來決定要不要買日馳。說不定明天就崩到底給你看XDD。

LPPL

網站都已經寫的滿清楚的,我就不多加贅述,雖然你第一次看可能會不太懂,可以先看以下的補充說明:一些關於model(模型)的觀念:

model就是一個數學模型,由一或多個數學公式組成,可以拿來預測天氣、預測人口成長率、預測各種東西,是個很強大的東西!

LPPL就是一個數學震盪公式:

y(t)=LPPL(t,tcrash,x)y(t) = LPPL(t, t^\text{crash}, x)

當中的一些變數:

  1. y(t)y(t):模擬的股價 1. tt:就是指每天的日期
  2. tcrasht^\text{crash}:崩盤的日期 3. xx:其它一些比較不重要的參數

我們的目標就是調整 tcrasht^\text{crash}xx 的值,讓模擬的股價 y(t)y(t) 可以對應到某支股票的實際股價 ym(t)y^m(t),就可以了。通常我們會用以下的方法對應到一個股價: mintcrash,xy(t)ym(t)2\min_{ t^\text{crash}, x}\left|y(t)-y^{m}(t)\right|^{2}

這個步驟就是調整參數,讓股價模型預測的股價趨近一樣,通常也可以稱為最佳化!
於是,我們就有了 tcrasht^\text{crash} 的值,也就可以預測崩盤的時機了!

接下來可以參考網站中的介紹,看一下LPPL的預測能力。

以下示範日馳給大家看:

取得日馳股價

首先我們用 pandas 中的 datareader,取得日馳的股價

取得股價
1
2
3
4
5
6
7
8
from pandas_datareader import data # pip install pandas_datareader
import matplotlib.pyplot as plt # pip install matplotlib
import pandas as pd # pip install pandas
%matplotlib inline

data = data.DataReader("1526.tw", "yahoo", "2017-04-01","2018-01-10")
c = data['Close']
c.plot()

最佳化

接下來,將網站中的LPPL複製過來,並且將日馳的股價 data.Close assign給 sz 這個變數。

這段code是在計算最佳參數,我沒有仔細閱讀code。有時寫程式就要學會複製貼上,可以用,就好了,對我們來說這不是科學問題,而是工程問題XD。

LPPL模型分析
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
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fmin_tnc
import random
import pandas as pd
from pandas import Series, DataFrame
import datetime
import itertools

sz = data.Close
date = sz.index
time = np.linspace(0, len(sz)-1, len(sz))
close = np.array(np.log(sz))
DataSeries = [time, close]

def lppl (t,x): #return fitting result using LPPL parameters
a = x[0]
b = x[1]
tc = x[2]
m = x[3]
c = x[4]
w = x[5]
phi = x[6]
return a + (b*np.power(tc - t, m))*(1 + (c*np.cos((w *np.log(tc-t))+phi)))

def func(x):
delta = [lppl(t,x) for t in DataSeries[0]]
delta = np.subtract(delta, DataSeries[1])
delta = np.power(delta, 2)
return np.sum(delta)

class Individual:
"""base class for individuals"""

def __init__ (self, InitValues):
self.fit = 0
self.cof = InitValues

def fitness(self):
try:
cofs, nfeval, rc = fmin_tnc(func, self.cof, fprime=None,approx_grad=True, messages=0)
self.fit = func(cofs)
self.cof = cofs
except:
#does not converge
return False

def mate(self, partner):
reply = []
for i in range(0, len(self.cof)):
if (random.randint(0,1) == 1):
reply.append(self.cof[i])
else:
reply.append(partner.cof[i])
return Individual(reply)

def mutate(self):
for i in range(0, len(self.cof)-1):
if (random.randint(0,len(self.cof)) <= 2):
self.cof[i] += random.choice([-1,1]) * .05 * i

def PrintIndividual(self):
#t, a, b, tc, m, c, w, phi
cofs = "A: " + str(round(self.cof[0], 3))
cofs += "B: " + str(round(self.cof[1],3))
cofs += "Critical Time: " + str(round(self.cof[2], 3))
cofs += "m: " + str(round(self.cof[3], 3))
cofs += "c: " + str(round(self.cof[4], 3))
cofs += "omega: " + str(round(self.cof[5], 3))
cofs += "phi: " + str(round(self.cof[6], 3))
return "fitness: " + str(self.fit) +"\n" + cofs

def getDataSeries(self):
return DataSeries

def getExpData(self):
return [lppl(t,self.cof) for t in DataSeries[0]]

def getTradeDate(self):
return date

def fitFunc(t, a, b, tc, m, c, w, phi):
return a - (b*np.power(tc - t, m))*(1 + (c*np.cos((w *np.log(tc-t))+phi)))

class Population:
"""base class for a population"""

LOOP_MAX = 1000

def __init__ (self, limits, size, eliminate, mate, probmutate, vsize):
'seeds the population'
'limits is a tuple holding the lower and upper limits of the cofs'
'size is the size of the seed population'
self.populous = []
self.eliminate = eliminate
self.size = size
self.mate = mate
self.probmutate = probmutate
self.fitness = []
for i in range(size):
SeedCofs = [random.uniform(a[0], a[1]) for a in limits]
self.populous.append(Individual(SeedCofs))

def PopulationPrint(self):
for x in self.populous:
print(x.cof)

def SetFitness(self):
self.fitness = [x.fit for x in self.populous]

def FitnessStats(self):
#returns an array with high, low, mean
return [np.amax(self.fitness), np.amin(self.fitness), np.mean(self.fitness)]

def Fitness(self):
counter = 0
false = 0
for individual in list(self.populous):
print('Fitness Evaluating: ' + str(counter) + " of " + str(len(self.populous)) + " \r"),
state = individual.fitness()
counter += 1
if ((state == False)):
false += 1
self.populous.remove(individual)
self.SetFitness()
print("\n fitness out size: " + str(len(self.populous)) + " " + str(false))

def Eliminate(self):
a = len(self.populous)
self.populous.sort(key=lambda ind: ind.fit)
while (len(self.populous) > self.size * self.eliminate):
self.populous.pop()
print("Eliminate: " + str(a- len(self.populous)))

def Mate(self):
counter = 0
while (len(self.populous) <= self.mate * self.size):
counter += 1
i = self.populous[random.randint(0, len(self.populous)-1)]
j = self.populous[random.randint(0, len(self.populous)-1)]
diff = abs(i.fit-j.fit)
if (diff < random.uniform(np.amin(self.fitness), np.amax(self.fitness) - np.amin(self.fitness))):
self.populous.append(i.mate(j))
if (counter > Population.LOOP_MAX):
print("loop broken: mate")
while (len(self.populous) <= self.mate * self.size):
i = self.populous[random.randint(0, len(self.populous)-1)]
j = self.populous[random.randint(0, len(self.populous)-1)]
self.populous.append(i.mate(j))
print("Mate Loop complete: " + str(counter))

def Mutate(self):
counter = 0
for ind in self.populous:
if (random.uniform(0, 1) < self.probmutate):
ind.mutate()
ind.fitness()
counter +=1
print("Mutate: " + str(counter))
self.SetFitness()

def BestSolutions(self, num):
reply = []
self.populous.sort(key=lambda ind: ind.fit)
for i in range(num):
reply.append(self.populous[i])
return reply;

random.seed()

from matplotlib import pyplot as plt
import datetime
import numpy as np
import pandas as pd
import seaborn as sns
sns.set_style('white')

limits = ([8.4, 8.8], [-1, -0.1], [350, 400], [.1,.9], [-1,1], [12,18], [0, 2*np.pi])
x = Population(limits, 20, 0.3, 1.5, .05, 4)
for i in range (2):
x.Fitness()
x.Eliminate()
x.Mate()
x.Mutate()

x.Fitness()
values = x.BestSolutions(3)
for x in values:
print (x.PrintIndividual())

別懷疑,就給它跑下去就對了,我也沒仔細看就跑了,需要跑一段時間。
跑完後應該會得到下圖:

可以看到這段代碼總共做了3次預測,以防算出來不準,其中的 critical time 就是指崩盤日期,約在380天左右,而我們總共測了200天,所以距離今天(2018-01-10),大概還需要 380 - 200 = 180 個交易日,才會有毀滅性的崩盤。

最後畫出圖,我也是複製貼上:

畫圖
1
2
3
data = pd.DataFrame({'Date':values[0].getDataSeries()[0],'Index':values[0].getDataSeries()[1],'Fit1':values[0].getExpData(),'Fit2':values[1].getExpData(),'Fit3':values[2].getExpData()})
data = data.set_index('Date')
data.plot(figsize=(14,8))

結論

LPPL這個模型看看就好,財報、月報也很重要,依照這個程式的結果,日馳距離泡沫崩盤還要很久,但千萬別因為這樣就買入,畢竟股票有太多要考量的因素了,大戶的脾氣不是數學模型可以模擬的XDDDDDDD…。

假如覺得文章不錯,那更不能錯過我們的影音課程喔!
或我們按個 鼓勵一下吧!