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

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

預防針

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

LPPL

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

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

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

當中的一些變數:

  1. y(t):模擬的股價 1. t:就是指每天的日期
  2. 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…。

FinLab - 韓承佑

嗨大家好,我是韓承佑,FinLab創辦人,畢業於巴黎薩克雷大學資工博士,目前擔任臺灣量化交易協會 學術顧問、台北商業大學 創新育成中心 創業技術顧問與上市科技公司 量化交易顧問。當初,我喜歡寫程式、無意間因為軟體比賽接觸Fintech,從此開始了財經跟程式的學習之路。我們成立 FinLab 量化投資部落格,用自己研發的軟體,對台灣股市做大量快速的實驗。希望可以在量化投資的路上,當大家的「武器製造商」!