# 破解股市泡沫之謎——對數周期冪率（LPPL）模型 > 來源：https://uqer.io/community/share/567a4fbd228e5b344568810f ## 引言 雖然離開物理專業有好幾年了，但一直有些念念不忘，碼著代碼，寫著開題報告，又閑不住想來講一個沒有好奇心的物理學家不是好金融學家的故事。 發現金融泡沫并預測到其何時破裂是很多從事金融行業的人的夢想。如今中國股市也成為了熱門的話題，然而，資本狂歡之后是股災，多少人因此從千萬富翁炒股變成百萬富翁，預測泡沫是所有人的夢想。 我們的主角Sornette教授登場了。Didier Sornette是一位受過培訓的統計物理學家和地球物理學家，目前在瑞士聯邦理工學院蘇黎世分校(Swiss Federal Institute of Technology in Zurich)任金融學教授，主講創業風險。他似乎并沒有因外界對這種綜合學科研究方法的熱情有所減弱而感到煩惱。相反，他還在做自己大部分職業生涯一直在做的事：不僅在主要物理期刊上發表文章，還在領先的金融期刊上發表文章。 Sornette教授開始嘗試著解答這個問題——不是通過傳統的金融學方法，而是將物理學思想引入其中。作為2004年出版的《股市為什么會崩盤》(Why Stock Markets Crash)一書的作者，Sornette教授實質上是希望更深刻地理解泡沫的形成和發展。在對復雜體系的分析中，他獨自——或者是和極少數幾個人一起——引領著三個并行領域：純物理學、應用經濟學和計量經濟學，以及市場從業人員。 在《股市為什么會崩盤》一書中，Sornette教授全面分析了一個由其提出的預測市場泡沫的模型——對數周期冪律（LPPL）模型。該模型對之后許多次市場泡沫都進行了準確的預測，由于該模型由Johansen，Ledoit和Sornette共同提出并完善，因此也被稱為JLS模型。我們來聊一聊它。 ## 什么是對數周期冪率模型？ 作為純物理學家的Sornettee教授不甘于僅僅在物理學領域有所建樹，他還看到了金光閃閃的華爾街，在那里，各類煉金術師在尋找各種允許少數人持續獲利的方法。于是，Sornettee教授在金融領域的跨界之旅開始了。他腦洞大開，想將物理學模型延伸到金融學領域中，而他找到的第一把金光閃閃的鑰匙叫做易辛模型——一種描述物質鐵磁性的經典模型。簡單地說，易辛模型認為單個原子的磁矩只可能有兩種狀態，+1（自旋向上）或者-1（自旋向下），原子以某種規則排列著，并存在著交互作用，使得相鄰之間的原子的自旋互相影響。 Sornette教授的眼睛仿佛一下子充滿了光芒，他仿佛看到了美元紙幣上的林肯在向他招手。受易辛模型啟發，Sornette教授認為在金融市場中，投資者也只具有兩種狀態，即買或者賣。同時，投資者的交易行為取決于其他投資者的決策及外部因素的影響，這與易辛模型是多么的相似！ 假想我們處于這樣的一個市場中：資產沒有派息、銀行利率為零、市場極度厭惡風險，并且市場有著充足的流動性。顯然，在這個市場中的金融資產沒有任何價值，也就是其基礎價值為零。在這樣的框架內，市場中出現兩類投資者，如上文所說，一類是理性投資者，一類是非理性的噪聲投資者。后者具有羊群效應，使得金融資產價格偏離其基礎價值，在沒有足夠的做空機制下，該結果導致理性投資者也不得不跟隨噪聲投資者的行為，通過享受泡沫來獲得收益。最終當趨勢達到某一臨界值時，大量投資者沒有足夠的頭寸維持該趨勢，于是手中的賣單導致了市場的崩盤。 那么這是一個怎樣的趨勢呢？Sornette教授考慮了自激勵的正反饋過程的思想，而該過程會導致大量交易者的行為方式逐漸趨于一致。在經過一些推導之后，Sornette教授發現該趨勢是按對數周期冪律（LPPL）增長，這里給出唯一的也是最重要的公式。  這里不去探討該公式的具體意義，讓我們看一下它的樣子。  可以看到，隨著時間增長，資產價格有著近似指數增長的特點，但同時也伴隨著不斷的振蕩，隨著時間越來越接近臨界時間，振蕩的幅度逐漸減小，增長速度逐漸增大，進入超指數的增長狀態，最終市場在臨界時間點附近崩盤。通過該模型，人們可以提前獲知可能的臨界時間點來規避風險。該模型曾成功預測了2008年的石油泡沫，美國房地產泡沫，以及2009 年中國股市泡沫等。 然而，試圖根據泡沫跡象采取行動的交易員，現在或許會非常失望。正如Sornette教授自己承認的，其理論實際上旨在估計這種泡沫的存在時間，而過早退出市場將是個錯誤，很可能會損失大量資金，然而離開過晚就不是幾個人失去工作的事了。事實上，該模型并沒有考慮交易者以外的因素，比如政策層面或者市場情緒等因素，但將金融系統認為是一個復雜系統并加以研究的思想是深遠的。現在也有許多將語義情緒分析等類似機器學習的方法應用于模型中來對市場狀態進行分析。 毫無疑問，更多地了解泡沫的形成和發展，價值是無法估量的。關注經濟和金融以外的領域有助于拓展思路，但不要指望找到一個指導市場交易的萬能公式。 LPPL模型收到的批評與收到的贊揚一樣多，有不少人認為該模型沒有操作價值。如果你想了解更多關于它的信息，可以仔細詳讀[Everything You Always Wanted to Know about Log Periodic Power Laws for Bubble Modelling but Were Afraid to Ask](https://kr.mathworks.com/matlabcentral/answers/uploaded_files/12237/frolova-paper-15.pdf)。 最后，如果你有好奇心，那么你一定想知道LPPL模型的實戰結果到底如何。下面我用優礦再現了LPPL模型預測2015年夏天A股市場的泡沫。 P.S 這次股災雖然我逃頂了，也許我馬后炮了~ :) 使用`DataAPI.MktIdxdGet()`函數獲取上證指數2014年1月1日至2015年6月10日的指數信息（股災發生在約一星期后）。`lib`庫我也放在了文章最后，大家可以嘗試著使用。特別感謝jd8001，用于擬合的GA算法的核心代碼框架來自于他，我對參數設置做了細致的調整，以讓它更好的符合A股市場，最后我友好地加上了一些注釋。 ```py import lib.lppltool as lppltool 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 = lppltool.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() Fitness Evaluating: 0 of 20 Fitness Evaluating: 1 of 20 Fitness Evaluating: 2 of 20 Fitness Evaluating: 3 of 20 Fitness Evaluating: 4 of 20 Fitness Evaluating: 5 of 20 Fitness Evaluating: 6 of 20 Fitness Evaluating: 7 of 20 Fitness Evaluating: 8 of 20 Fitness Evaluating: 9 of 20 Fitness Evaluating: 10 of 20 Fitness Evaluating: 11 of 20 Fitness Evaluating: 12 of 20 Fitness Evaluating: 13 of 20 Fitness Evaluating: 14 of 20 Fitness Evaluating: 15 of 20 Fitness Evaluating: 16 of 20 Fitness Evaluating: 17 of 20 Fitness Evaluating: 18 of 20 Fitness Evaluating: 19 of 20 fitness out size: 20 0 Eliminate: 14 Mate Loop complete: 25 Mutate: 2 Fitness Evaluating: 0 of 31 Fitness Evaluating: 1 of 31 Fitness Evaluating: 2 of 31 Fitness Evaluating: 3 of 31 Fitness Evaluating: 4 of 31 Fitness Evaluating: 5 of 31 Fitness Evaluating: 6 of 31 Fitness Evaluating: 7 of 31 Fitness Evaluating: 8 of 31 Fitness Evaluating: 9 of 31 Fitness Evaluating: 10 of 31 Fitness Evaluating: 11 of 31 Fitness Evaluating: 12 of 31 Fitness Evaluating: 13 of 31 Fitness Evaluating: 14 of 31 Fitness Evaluating: 15 of 31 Fitness Evaluating: 16 of 31 Fitness Evaluating: 17 of 31 Fitness Evaluating: 18 of 31 Fitness Evaluating: 19 of 31 Fitness Evaluating: 20 of 31 Fitness Evaluating: 21 of 31 Fitness Evaluating: 22 of 31 Fitness Evaluating: 23 of 31 Fitness Evaluating: 24 of 31 Fitness Evaluating: 25 of 31 Fitness Evaluating: 26 of 31 Fitness Evaluating: 27 of 31 Fitness Evaluating: 28 of 31 Fitness Evaluating: 29 of 31 Fitness Evaluating: 30 of 31 fitness out size: 31 0 Eliminate: 25 Mate Loop complete: 25 Mutate: 0 Fitness Evaluating: 0 of 31 Fitness Evaluating: 1 of 31 Fitness Evaluating: 2 of 31 Fitness Evaluating: 3 of 31 Fitness Evaluating: 4 of 31 Fitness Evaluating: 5 of 31 Fitness Evaluating: 6 of 31 Fitness Evaluating: 7 of 31 Fitness Evaluating: 8 of 31 Fitness Evaluating: 9 of 31 Fitness Evaluating: 10 of 31 Fitness Evaluating: 11 of 31 Fitness Evaluating: 12 of 31 Fitness Evaluating: 13 of 31 Fitness Evaluating: 14 of 31 Fitness Evaluating: 15 of 31 Fitness Evaluating: 16 of 31 Fitness Evaluating: 17 of 31 Fitness Evaluating: 18 of 31 Fitness Evaluating: 19 of 31 Fitness Evaluating: 20 of 31 Fitness Evaluating: 21 of 31 Fitness Evaluating: 22 of 31 Fitness Evaluating: 23 of 31 Fitness Evaluating: 24 of 31 Fitness Evaluating: 25 of 31 Fitness Evaluating: 26 of 31 Fitness Evaluating: 27 of 31 Fitness Evaluating: 28 of 31 Fitness Evaluating: 29 of 31 Fitness Evaluating: 30 of 31 fitness out size: 31 0 fitness: 0.99612688166 A: 9.817B: -0.681Critical Time: 365.323m: 0.207c: -0.023omega: 12.241phi: 4.25 fitness: 0.99612688166 A: 9.817B: -0.681Critical Time: 365.323m: 0.207c: -0.023omega: 12.241phi: 4.25 fitness: 0.99653502204 A: 9.8B: -0.667Critical Time: 365.405m: 0.209c: -0.023omega: 12.267phi: 4.105 ``` ```py 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)) <matplotlib.axes.AxesSubplot at 0x663c250> ```  模型預測的臨界時間（Critical Time）為365，即為6月10日（350）之后的15個交易日左右。實際股災時間為6月15日（353），比實際結果晚10個交易日左右。 `lib`庫代碼，請保存并命名為`lppltool` ```py #code created by jd8001 #reference: https://github.com/jd8001/LPPL #kindly thank jd8001! 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 SP = DataAPI.MktIdxdGet(ticker='000001',beginDate='20140101',endDate='20150610',field=["tradeDate","closeIndex"],pandas="1") global date = SP.tradeDate time = np.linspace(0, len(SP)-1, len(SP)) close = [np.log(SP.closeIndex[i]) for i in range(len(SP))] global DataSeries 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]] #生成lppl時間序列 delta = np.subtract(delta, DataSeries[1]) #將生成的lppl時間序列減去對數指數序列 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): # 交配，0.5的概率自身的參數保留，0.5的概率留下partner的參數，即基因交換 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): #print "Mutate" + str(i) 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 #return str(self.cof) + " fitness: " + str(self.fit) 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() ```