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

預防針
這邊還是要提醒大家,這個只是個有趣的小實驗,但我真的不確定這個很可靠,請不要依賴LPPL買賣股票,可以當作一個參考,追蹤一下日馳的股價是否跟LPPL模型顯示的一樣。我個人是持非常懷疑的態度!請大家不要因為這篇文章來決定要不要買日馳。說不定明天就崩到底給你看XDD。
LPPL
該網站都已經寫的滿清楚的,我就不多加贅述,雖然你第一次看可能會不太懂,可以先看以下的補充說明:一些關於model(模型)的觀念:
model就是一個數學模型,由一或多個數學公式組成,可以拿來預測天氣、預測人口成長率、預測各種東西,是個很強大的東西!
LPPL就是一個數學震盪公式:

當中的一些變數:
- y(t):模擬的股價 1. t:就是指每天的日期
 - tcrash:崩盤的日期 3. x:其它一些比較不重要的參數
 
我們的目標就是調整 tcrash 和 x 的值,讓模擬的股價 y(t) 可以對應到某支股票的實際股價 ym(t),就可以了。通常我們會用以下的方法對應到一個股價:

這個步驟就是調整參數,讓股價跟模型預測的股價趨近一樣,通常也可以稱為最佳化!
於是,我們就有了 tcrash 的值,也就可以預測崩盤的時機了!
接下來可以參考網站中的介紹,看一下LPPL的預測能力。
以下示範日馳給大家看:
取得日馳股價
首先我們用 pandas 中的 datareader,取得日馳的股價
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模型分析
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 個交易日,才會有毀滅性的崩盤。
最後畫出圖,我也是複製貼上:
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…。