来源：华泰金融工程
林晓明 S0570516010001 研究员
陈 烨 S0570518080004 研究员
李子钰 S0570519110003 研究员
何 康 S0570118080081 联系人
源洁莹 S0570119080125 联系人
王晨宇 S0570119110038 联系人
报告发布时间：2020年5月28日
摘要
藏本模型是刻画耦合振子同步现象的经典模型
藏本模型是刻画耦合振子同步现象的经典模型。同步现象在自然界和人类社会广泛存在，其本质是系统中的微观个体相互耦合最终形成宏观级别的周期。藏本模型的基本形式是描述系统中每个振子相位的一维动力学方程。模型中，振子受到其它振子的吸引而趋于同步，其等价形式是振子受到平均场的吸引而趋于同步。振子数量N、耦合强度K、自然频率洛伦兹分布半高宽γ这三项关键参数对同步的影响不尽相同。Sakaguchi－Kuramoto模型是藏本模型的变式，可借助该模型研究振子合作、叛逆等行为对系统终态的影响。
同步现象无处不在，本质是系统中的微观个体相互耦合形成宏观周期
同步现象无处不在，心脏起搏细胞自发同步产生有节律的心跳，夏夜萤火虫的同步闪烁，演出结束后观众自发产生有节奏的掌声，原子共同受激发射产生激光，钟摆或节拍器的同步摆动，全球经济指标及金融资产的同起同落，这些看似神奇的现象背后存在共通之处。系统的每个微观个体可视作周期运动的振子，众多个体通过相互间的耦合，最终形成步调一致的宏观周期，物理学中称为耦合振子同步。藏本模型是刻画耦合振子同步现象的经典模型。
藏本模型的基本形式是关于振子相位的一维动力学方程
藏本模型的基本形式是描述系统中每个振子相位的一维动力学方程。振子相位的变化量（即运动速度）由1）该振子本身自然频率ω，2）耦合强度K，3）该振子和其它振子的相位差决定。耦合强度K＞0使得振子间相互吸引。引入序参量、平均相位和平均场的概念后，振子受到其它振子的吸引而趋于同步，等价于振子受到平均场的吸引而趋于同步。
藏本模型中三项关键参数对同步的影响不尽相同
藏本模型中振子数量N、耦合强度K、自然频率洛伦兹分布半高宽γ这三项关键参数对同步的影响不尽相同。当振子数量N≥3时，振子数量对同步速度影响不大。耦合强度K小于某个定值时，无法同步；大于某个定值时，K越大同步越快。振子自然频率服从洛伦兹分布的半高宽γ越大，系统越难同步。当振子的自然频率在某个特定范围内时，振子能与平均场同步；当振子的自然频率超出该范围时，振子无法同步。
借助藏本模型变式研究振子合作、叛逆等行为对系统终态的影响
Sakaguchi－Kuramoto模型是藏本模型的变式，可借助该模型研究振子合作、叛逆等行为对系统终态的影响。当系统由合作和叛逆振子组成时：若合作振子占优，则两类振子分别形成团簇而同步，团簇间相位差为π，即处于反相位；若叛逆振子占优，则系统无法同步。当系统由合作和孤立振子组成时，孤立振子始终按自身相位运动，合作振子形成同步。当系统由合作和延迟振子组成时，两类振子分别形成团簇而同步，团簇间存在相位延迟。合作振子占比越高，系统越倾向于同步；根据不同的合作振子占比p值，存在非相干态、π态和行波态三种不同的系统状态。
风险提示：包括藏本模型在内的动力学系统模型是对真实世界的简单刻画，存在过度简化的可能；周期规律基于历史数据总结，历史规律可能失效；周期规律对市场长期规律进行研究，无法判断短期的市场情绪与政策冲击带来的影响；市场在某些极端情形下可能出现规律以外的交易风险。
藏本模型及其变式
耦合振子同步现象
同步现象无处不在。人类生命的维系离不开心脏近一万个起搏细胞的同步振荡，如果各细胞放电节律不一致，就可能出现心律不齐。夏夜林中萤火虫的同步闪烁，释放出强烈的求偶信号，使整个种群得以延续。演出大厅内，观众自发产生有节奏的鼓掌声。全球经济指标及金融资产时常同起同落。即便在非生命的世界，原子共同受激发射产生激光，钟摆或节拍器的同步摆动等等，这些同步现象也早已为人们熟知。
上述现象背后的共同点是，系统的每个微观个体都可视作周期运动的振子，众多个体通过相互间的耦合，最终形成步调一致的宏观周期，在物理学中称为耦合振子同步。奇妙的现象引发人们思考：同步背后的原理是什么？有没有模型能够刻画这些现象？
在华泰金工《周期起源2：周期趋同现象的动力学系统模型》（20200102）一文中，我们展示了如何通过Haken-Kelso-Bunz模型（简称HKB模型）解释心理学中的手指运动周期趋同现象，并提出我们对经济系统周期现象起源的一种猜测。在复杂系统领域，还有一些模型也能刻画耦合振子同步现象，应用范围较HKB模型更广，其中最具代表性的模型之一是藏本模型（Kuramoto Model）。
本文我们将首先介绍藏本模型及其变式Sakaguchi－Kuramoto模型和Hong－Strogatz模型（其中Hong－Strogatz模型可视作Sakaguchi－Kuramoto模型的特例）。随后我们将对藏本模型和Sakaguchi－Kuramoto模型进行模拟，展示耦合振子同步现象，并探讨关键参数以及振子合作、叛逆等行为对同步结果的影响。
基础知识
振子在空间内按一定规律振动，如常见的简谐振动，其运动轨迹可以由直角坐标系下的正弦波刻画，如左下图所示。以极坐标系视角考察，如右下图所示，振子在圆周上沿逆时针方向做周期运动，振子所处方向和极轴之间的夹角为相位θ，单位为弧度（rad）；运动的速度以角频率ω表示，即单位时间内相位的变化量θ ̇，等价于相位的一阶导数，也称为振子的固有频率或自然频率，单位为弧度/秒（rad/s）。
藏本模型
藏本模型由日本物理学家藏本由纪（Yoshiki Kuramoto）于1975年提出，形式简洁且易于理解，是描述耦合振子同步现象的经典模型，广泛应用于各领域。藏本模型的基本形式是描述系统中每个振子相位的一维动力学方程：
其中，状态变量θi代表振子i的相位。等号左侧为θi的一阶导数，即相位的变化量，反映振子的运动速度。参数ωi代表振子i的固有频率，一般假定服从洛伦兹分布、高斯分布等，概率密度函数以g(ω)表示。参数N为系统中振子总数。参数K代表耦合强度，即振子之间的相互作用强度，一般取K＞0，即振子之间相互“吸引”。求和符号Σ内的sin(θj－θi)代表振子j对振子i的作用，整个Σ代表系统内全部振子对振子i的影响之和。
藏本模型通俗化理解的经典案例是演出大厅里观众鼓掌声的同步现象，在Neda等（2000）和鞠萍（2016）文献中均有介绍。演出结束时，剧场内N个观众开始鼓掌，每个人自身的鼓掌速度ωi有快有慢不尽相同，每个人鼓掌的相位θi也不同，在宏观上表现为杂乱无章的掌声。理想状态下，每个人会受到其他人掌声的“吸引”，掌声相位的变化量受相位差θj－θi影响，相位θi不断作出调整。随着时间的进行，掌声最终趋于同步，每个人的相位θi保持一致。
下面我们进一步阐述振子之间的相互作用机制。假设振子i和振子j的相位差θj－θi介于-π到π之间（若超过此范围，相当于长跑比赛中的“套圈”，可通过加减2π的整数倍而转化到此范围）：
1.当θj－θi＞0，即振子i落后振子j时，如左上图所示，sin(θj－θi)＞0，由于K＞0，此时振子i相位的变化量等于其固有频率加上一个正数，那么下一时刻振子i的相位将相对增大，即“加速追赶”振子j。
2.当θj－θi＜0，即振子i领先振子j时，如右上图所示，sin(θj－θi)＜0，由于K＞0，此时振子i相位的变化量等于其固有频率加上一个负数，那么下一时刻振子i的相位将相对减小，即“减速等待”振子j。
3.当θj－θi＝0，即振子i和振子j同步时，sin(θj－θi)＝0，此时振子i相位的变化量等于其固有频率，振子j对振子i无影响，是一种相对稳定的状态。
总的来看，K＞0使得振子之间存在相互“吸引”作用。直观地看，吸引作用才能使得各振子最终趋于同步。
Sakaguchi－Kuramoto模型
藏本模型较为简洁，后续的研究者提出各种藏本模型的变式，以适应更复杂的研究场景。1986年Sakaguchi和Kuramoto提出Sakaguchi－Kuramoto模型，引入相位延迟β，此时动力学方程为：
其中，除β外的各个变量含义均与原版藏本模型相同。相位延迟β的含义是，振子i并非倾向于和其它任意振子j同步，而是保持β的相位差。当θi＝θj－β，即振子i落后振子j的相位等于β时，sin(θj－θi－β)＝0，此时振子j对振子i无影响，是一种相对稳定的状态。在K＞0的假设下，习惯上β的范围取－π/2至π/2之间。
Hong－Strogatz模型
2011年Hong和Strogatz提出Hong－Strogatz模型，该模型和原始藏本模型的不同之处在于：原始藏本模型中，每个振子受其它振子的影响程度相同，耦合强度均为K；Hong－Strogatz模型中，各振子受其它振子的影响程度不同，耦合强度以参数Ki表示：
原始藏本模型的K一般取正数，即认为振子之间为相互吸引关系。Hong－Strogatz模型中，Ki拓展至任意实数。当振子i的耦合强度Ki＞0时，该振子倾向于受其它振子吸引，相当于系统中的“合作者”；当振子i的耦合强度Ki＜0时，该振子倾向于和其它振子相排斥，相当于系统中的“叛逆者”。
特别地，当Ki＜0时，“叛逆”振子的动力学方程可以写作：
振子i倾向于和其它任意振子j保持π的相位差，即处于反相位的状态，正好相当于Sakaguchi－Kuramoto模型中相位延迟β取π的情形。因此，一般将Hong－Strogatz模型视作Sakaguchi－Kuramoto模型的特例。当Sakaguchi－Kuramoto模型允许“叛逆”振子的存在时，相位延迟β的范围可拓展到-π至π之间。
复序参量、序参量和平均相位
下面两节将介绍藏本模型及其变式中的重要概念。为了更清晰地描述振子的动力学行为，引入下式：
Z定义为系统内各振子eiθ的算数平均值，i代表虚单位。由上式可知Z为复数，可能包含虚部，称Z为耦合振子系统的复序参量（Complex Order Parameter）。对于复数Z，可以表示成reiψ的形式，其中r为复数Z的模（Modulus），也称为系统的序参量（Order Parameter）；ψ为复数Z的幅角（Argument），也称为系统的平均相位（Average Phase）。
下面我们从几何角度形象化地理解上述概念。首先需要明确，相位θ是针对单个振子而言，属于微观量；复序参量Z、序参量r和平均相位ψ是针对整个系统而言，属于宏观量。如上图所示，假设系统内各振子都位于单位圆上，每个振子都可以视作一个单位向量，将这些向量相加并求平均，即得到描述整个系统的向量Z，即复序参量。
Z相当于各个振子产生“合力”的均值，Z的长度对应模长r，即序参量；Z与极轴的夹角对应幅角ψ，即平均相位。模长r应小于1，由于各振子“发力”方向不同，与Z垂直方向的分量相互抵消，只保留Z方向的分量，因此“合力”小于每个分力的均值。当且仅当所有振子完全同步时，序参量r为1，此时所有振子“劲往一处使”。如果所有振子均匀分布在单位圆上，序参量r为0。一般而言，序参量r介于0和1之间，越接近1代表系统的同步效果越好。
原始藏本模型中，振子的相位改变量由1）该振子本身自然频率，2）振子耦合强度，3）该振子和其它所有振子的相位差决定。引入序参量r和平均相位ψ的概念后，可以将藏本模型进行恒等变换，得到如下形式：
此时，振子的相位改变量由1）该振子本身自然频率，2）振子耦合强度，3）耦合振子系统的序参量和平均相位决定。换言之，振子受到其它振子的吸引而趋于同步，等价于振子受到“平均场”的吸引而趋于同步。恒等变换的推导过程请参见附录。
平均场频率、旋转坐标系和有效频率
复序参量Z的本质是将各个振子取平均，用以描述整个系统的状态。由于各个振子在极坐标系下沿圆周做周期运动，那么其“均值”Z（即整个系统）也在做周期运动。定义Z的幅角ψ（即系统的平均相位）的一阶导数为系统的平均场频率Ω，即系统振动的角频率，反映系统的运动速度。
接下来我们将进行参照系的转换。不妨以体育场内的长跑比赛为例，说明参照系转换：
1. 以观众为参照系：所有运动员都在绕跑道向前运动，领先的运动员速度最快（角频率最高），落后的运动员速度最慢（角频率最低）；此时参照系为静止坐标系。
2. 以最后一位运动员为参照系：除该运动员以外的运动员都在绕跑道向前运动，速度为其它运动员和自身的速度差；该运动员自身处于静止状态；此时参照系处于运动状态，运动速度等于最后一位运动员的速度，可视作旋转坐标系。
3. 以所有运动员的平均运动速度为参照系：速度靠前的运动员在向前运动，速度靠后的运动员反而向后运动；注意此时参照系也处于运动状态，仍可视作旋转坐标系。
仿照上面的第3种情况，我们将耦合振子系统的静止坐标转换为旋转坐标系。如上图所示，以复序参量Z作为参照系：Z的绝对运动速度为平均场频率Ω，以自身作为参照系后，Z的相对运动速度为0；振子i的绝对运动速度为θ ̇_i，以Z作为参照系后，振子i的相对运动速度为θ ̇_i－Ω。
由上节知，静止坐标系下，原始藏本模型可以写成序参量r和平均相位ψ的形式：
其中，平均相位ψ以平均场频率Ω旋转。改为旋转坐标系后，原始藏本模型就可以写成如下形式，其中，相位θ和频率ω均为旋转坐标系下的取值（详细推导过程请参见附录）：
为什么要做参照系转换？我们希望更便捷地判断系统是否同步。当全部振子实现同步时，各振子的相位相同，并且相位的一阶导数也相同，都等于系统的平均场频率。如果以Z作为参照系，各振子的相对运动速度就等于0，即相位的一阶导数等于0。记振子i在一段时间内旋转坐标系下相位θi一阶导数的均值为该振子的有效频率ωeff：
其中，<>t表示在一段时间内取均值。若系统完全同步，则各振子的有效频率应为0。
总的来看，我们得到判断系统是否完全同步的三种方式：
1. 绘制各振子运动相位，观察相位是否重合。
2. 从宏观角度看，整个系统的序参量r是否等于1。
3. 从微观角度看，每个振子的有效频率ωeff是否等于0。
藏本模型模拟
本章我们将对藏本模型进行模拟，观察耦合振子同步现象，并考察模型各项参数对同步的影响。回顾上一章介绍的原始藏本模型，N个振子构成的耦合振子系统中，任意振子i的动力学方程为：
若振子i的自然频率ωi服从洛伦兹分布，概率密度函数g(ω)表示为：
振子数量N、耦合强度K、自然频率ω服从洛伦兹分布的半高宽γ是最重要的三项参数。我们分别测试各项参数对同步的影响，本章核心结论如下图。
本章模拟所采用的具体参数如下表。需要说明的是：
1.前2组测试假设各振子自然频率均为ω，即常量；第3组测试假设各振子自然频率服从半高宽为γ的洛伦兹分布，即随机变量。
2.各组测试中，假设振子的初始相位等距分布在0至2π的圆周上，第i－1个振子的初始相位为2iπ/N。
3.前2组模拟采用静止坐标系，模拟时间间隔Δt为0.01；第3组模拟采用旋转坐标系，模拟时间间隔Δt为0.1，根据最后2000步模拟的振子运动速度计算有效频率。
振子数量N的影响
本节考察振子数量N对同步的影响。这里坐标系采用静止坐标系。
双振子情形（N＝2）
首先考察最简单的情形，假设系统内仅有2个振子，初始相位分别为0和π，即位于圆周上对称的两点，耦合强度K为1，即相互吸引关系，各振子自然频率ω均为0.1。
左下图为两个振子相位θ随时间变化情况。模拟开始阶段，相位相差π，当时间在20~40之间时，两个振子相位开始趋同，最终相位保持完全一致，实现同步。
为了更形象地展示振子的运动情况，对相位取正弦值sin(θ)，右下图为正弦值随时间变化情况。由图可知，两个振子在开始阶段各完整地走完了半个正弦波，随后振子2受振子1影响拐头向下，实现同步。同步后整个系统的周期运动频率等于各振子的自然频率（ω＝0.1）。
双振子系统序参量如左下图所示，平均相位和平均相位正弦值如右下图所示。开始阶段，序参量r为0，系统处于不同步状态，随后迅速上升，达到1。系统序参量等于1是完全同步的标志。
五振子情形（N＝5）
下面考察五振子的情形。各振子相位正弦值如左下图所示，开始阶段振子按各自相位运动，在时间为60左右趋于同步，最终相位正弦值完全重合。系统序参量和平均相位正弦值如右下图所示，序参量r开始为0，在时间为60左右开始迅速上升，达到1，实现同步。
不同振子数量N比较
测试不同振子数量N对同步的影响。左下图展示N＝2、5、10、20四种情形的序参量随时间变化情况，除N以外的模拟参数均保持一致。由图可知，N＝2时同步相对最快，N＝5、10、20三条线几乎重合。右下图展示N＝2~20时系统同步所需时间，定义序参量r超过0.99即完成同步。由图可知，N＝2时同步时间在30~40之间，N≥3时同步时间均在65~71之间，变化不大。由此得出结论：当振子数量N≥3时，振子数量对同步速度影响不大。
耦合强度K的影响
本节考察耦合强度K对同步的影响。这里坐标系采用静止坐标系。
强耦合情形（K＝0.8）
首先考察强耦合的情形，假设系统包含5个振子，初始相位等距分布在圆周上，耦合强度K为0.8。系统各振子相位正弦值如左下图所示，序参量和平均相位正弦值如右下图所示。观察可知同步的发生时间点在80左右，稍慢于上节K＝1的情形。
中等耦合情形（K＝0.5）
其次考察中等耦合的情形，耦合强度K为0.5，其余参数保持不变。系统各振子相位正弦值如左下图所示，序参量和平均相位正弦值如右下图所示。观察可知同步的发生时间点在110~130区间内，慢于此前K＝0.8的情形。
弱耦合情形（K＝0.2）
再次考察弱耦合的情形，耦合强度K为0.2，其余参数保持不变。观察可知同步的发生时间点在270~310区间内，慢于此前K＝0.5的情形。
不耦合情形（K≤0）
当耦合强度K为0时，振子间不存在互相影响；当耦合强度K为负数时，振子间为相互排斥关系。左下图和右下图分别展示K为0和-1情形下各振子相位正弦值。此时振子按各自相位运动，不发生同步。
不同耦合强度K比较
测试不同耦合强度K对同步的影响。左下图展示K＝0.2、0.5、0.8、1.0四种情形的序参量随时间变化情况，除K以外的模拟参数均保持一致。由图可知，这四种情况下，耦合强度越大，同步速度越快。
右下图展示K＝0~1时系统同步所需时间。由图可知，当K＜0.15时，耦合强度过小，无法同步；当K≥0.15时，同步所需时间随K的增加而降低。由此得出结论：耦合强度K小于某个定值时，无法同步；大于某个定值时，K越大同步越快。该定值的取值和其它参数有关，我们将在本章的最后讨论。
自然频率ω服从洛伦兹分布半高宽γ的影响
本节考察自然频率ω服从洛伦兹分布半高宽γ对同步的影响。本节和前两节测试相比，关键的不同点在于，振子数量N从5扩展到100，并且采用旋转坐标系，即以整个系统的平均场作为参照系。
由于振子数量大幅提升，考虑到运算时间开销，我们将模拟时间间隔Δt从0.01改为0.1。同样由于振子数量大幅提升，观察每个振子相位正弦值随时间变化以确定是否同步较困难，我们计算每个振子最后2000步模拟的相对运动速度均值作为有效频率。若有效频率等于0，说明振子相对于平均场的运动速度为0，即振子和平均场同步。
洛伦兹分布也称为柯西分布，是一种典型的肥尾分布。其概率密度函数可以写作：
不同半高宽γ下的洛伦兹分布概率密度函数图像如下所示。本节测试中，振子的自然频率将从该分布随机采样。为了避免极端值的影响，我们采用截断洛伦兹分布，确保采样得到的ω在±2之间。也有文献采用高斯分布，分布的选择不影响最终结论。
γ＝0.05情形
当γ为0.05，即自然频率取自较窄的分布时，各振子的自然频率接近0。系统序参量和平均相位正弦值随时间变化情况如左下图所示，这里展示时间0~100区间内的结果。由图可知，序参量接近1，平均相位接近0，说明系统接近完全同步。但由于少数离群振子的存在，系统未能达到完全同步。
系统内100个振子自然频率和有效频率如右上图所示。图中每个点代表一个振子，横轴为自然频率，纵轴为有效频率，即最后2000步模拟的平均速度。自然频率是振子固有的运动速度，相当于振子的“本性”；有效频率是振子在系统中受平均场持续影响，最终表现出的运动速度，相当于振子在集体中的行为。
由右上图可知，100个振子中99个振子的有效频率为0，即和平均场同步。但存在1个离群振子，该振子的自然频率为1.22，其绝对值高于其它全部振子，其有效频率为0.96。换言之，该振子本性“桀骜不驯”，最终未被平均场同化，没有和其它振子同步，但仍然一定程度上受到平均场的影响。
γ＝0.1情形
当γ为0.1时，系统序参量和平均相位正弦值随时间变化情况如左下图所示，系统内100个振子自然频率和有效频率如右下图所示。序参量r接近0.9，平均相位正弦值仍接近0，说明大多数振子实现同步。但右下图可以观察到7个振子的有效频率不为0，离群振子数量较γ＝0.05时更多。
我们推测，当振子自然频率超过一定临界值时，该振子不再和平均场同步。当γ为0.1时，能够同步的振子最大自然频率绝对值为0.68，无法同步的振子最小自然频率绝对值为0.88，推测此时临界值的绝对值在0.68和0.88之间。
γ＝0.2情形
当γ为0.2时，系统序参量和平均相位正弦值随时间变化情况如左下图所示，系统内100个振子自然频率和有效频率如右下图所示。序参量r在0.7附近波动，平均相位正弦值接近0，说明多数振子实现同步。但右下图可以观察到不少振子的有效频率不为0，离群振子数量较γ＝0.1时更多。
当γ为0.2时，能够同步的振子最大自然频率绝对值为0.59，无法同步的振子最小自然频率绝对值为0.63，推测此时振子同步自然频率临界值的绝对值在0.59和0.63之间。
γ＝0.5情形
当γ为0.5，即自然频率取自较宽的分布时，系统序参量和平均相位正弦值随时间变化情况如左下图所示，系统内100个振子自然频率和有效频率如右下图所示。序参量r接近0，平均相位正弦值波动较大，说明大多数振子不同步。右下图可以观察到仅有少数振子的有效频率为0。
当γ为0.5时，能够同步的振子最大自然频率绝对值为0.11，无法同步的振子最小自然频率绝对值为0.12，推测此时振子同步自然频率临界值的绝对值在0.11和0.12之间。
不同自然频率ω服从洛伦兹分布半高宽γ比较
测试不同自然频率ω服从洛伦兹分布半高宽γ对同步的影响。下图展示γ＝0.05、0.1、0.2、0.5四种情形的序参量随时间变化情况，除γ以外的模拟参数均保持一致。由图可知，γ越大，系统序参量越小，越难同步。
左下图展示γ＝0.05~0.5时，系统最后2000步模拟的序参量均值。γ越大，系统序参量越小，越难同步。右下图展示γ＝0.05~0.5时，系统内各振子有效频率绝对值的均值。γ越大，振子有效频率越大，越难同步。总的来看，振子自然频率服从洛伦兹分布的半高宽γ越大，系统越难同步。
临界频率的确定
通过以上测试，我们发现振子能否与平均场同步，取决于其自然频率，若自然频率在一定范围内，则能够同步，否则无法同步。而该临界频率又与系统有关，确切地说与全部振子自然频率服从分布的参数有关，当洛伦兹分布的半高宽γ较小时，临界频率较大。那么，能否定量计算该临界频率？
实际上，上述问题的答案就在藏本模型的动力学方程中。旋转坐标系下，藏本模型可以写成如下形式：
若振子和平均场同步，则振子相对于平均场的运动速度为0，即相对于平均场的相位θ的变化量为0，即等号左侧为0：
那么等号右侧也应为0，即：
若振子i的自然频率满足|ω_i |≤Kr，由于sin(θ_i)的取值范围为-1到1，一定存在某个θ_i的取值使得ω_i=Kr sin(θ_i)，此时θ ̇_i=0，振子实现同步。
若振子i的自然频率满足|ω_i |>Kr，由于sin(θ_i)的取值范围为-1到1，不存在任何θ_i的取值使得ω_i=Kr sin(θ_i)，此时θ ̇_i≠0，振子无法同步。
由此可得，临界频率应为耦合强度K和系统序参量r的乘积。当振子的自然频率介于-Kr和Kr之间时，振子能够和平均场同步；当振子的自然频率超出该范围时，振子无法同步。我们将临界频率的理论值Kr和本节实验观察得到的临界频率（即同步振子自然频率绝对值最大值max|ω|）相比较，如下表所示，最右侧两列基本匹配。
小结
本章探讨原始藏本模型参数对同步的影响，核心结论如下：
1. 当振子数量N≥3时，振子数量对同步速度影响不大。
2. 耦合强度K小于某个定值时，无法同步；大于某个定值时，K越大同步越快。
3. 振子自然频率服从洛伦兹分布的半高宽γ越大，系统越难同步。
4. 当振子的自然频率介于-Kr和Kr之间时，振子能够和平均场同步；当振子的自然频率超出该范围时，振子无法同步。
Sakaguchi－Kuramoto模型模拟
原始藏本模型包含两项重要假设：1）所有振子不存在相位延迟；2）所有振子耦合强度一致。Sakaguchi－Kuramoto模型和Hong－Strogatz模型分别针对以上两点做出改进，引入相位延迟β以及每个振子的耦合强度Ki。当振子的耦合强度为负时，相当于引入相位延迟β＝π。因此，Hong－Strogatz模型可以视作Sakaguchi－Kuramoto模型的特例。
本章我们将对Sakaguchi－Kuramoto模型进行模拟，关注以下问题：当系统中存在两类振子，一类为数量占比p的合作振子（受平均场吸引），另一类为数量占比（1－p）的排斥/孤立/延迟振子时，系统将呈现何种状态？合作振子的占比p如何影响系统最终状态？本章核心结论如下图所示。
本章模拟所采用的具体参数如下表所示。
两类振子组成的Sakaguchi－Kuramoto模型
本节测试两类振子组成的Sakaguchi－Kuramoto模型。假设系统包含5个振子，其中第一类振子为合作振子，耦合强度K为1，行为表现为受平均场吸引，数量占比为p；第二类振子耦合强度K＜1，数量占比为1－p。第二类振子的耦合强度K取不同值时，可能表现出叛逆、孤立、延迟等不同行为。
合作振子和叛逆振子
当第二类振子耦合强度K＜0时，振子行为表现为与平均场相排斥，不妨称为叛逆振子。首先我们测试合作振子和叛逆振子组成的系统，当合作振子数量占优时的系统状态。假设第一类合作振子的数量占比p为0.6，即包含3个振子，第二类叛逆振子耦合强度K为-1，数量占比为0.4，即包含2个振子。
各振子相位正弦值随时间变化情况如下图所示。开始阶段，各振子按各自的相位做周期运动，随后振子间的相互吸引或排斥作用逐渐体现，在时间220~240区间内发生了关键转变，最终3个第一类合作振子聚集成团簇（Cluster）从而实现同步，2个第二类叛逆振子聚集成另一个团簇从而也实现同步，两个团簇间的相位差为π，即相位相反。换言之，系统最终状态为“物以类聚，人以群分”，合作振子和叛逆振子分别形成各自的“阵营”。
为了更清晰地观察同步发生的细节，我们展示时间200~250区间内各振子的相位正弦值，如下图所示。可知3个第一类合作振子形成同步，2个第二类叛逆振子形成同步。
下面我们测试合作振子和叛逆振子组成的系统，当叛逆振子数量占优时的系统状态。假设第一类合作振子的数量占比p为0.4，即包含2个振子，第二类叛逆振子耦合强度K为-1，数量占比为0.6，即包含3个振子。
各振子相位正弦值随时间变化情况如下图所示。观察到各振子始终按各自的相位做周期运动，无法实现同步。当系统内叛逆振子占优时，即使是合作振子内部也无法聚集成团簇而实现同步。通俗地理解，当社会中破坏团结的个体过多，那么维护团结的个体也很难形成合力。
合作振子和孤立振子
下面我们测试合作振子和孤立振子组成的系统。孤立振子的耦合强度K为0，既不受平均场吸引，也不与平均场排斥。
首先测试合作振子占优的情形。假设第一类合作振子数量占比p为0.6，即包含3个振子，第二类孤立振子耦合强度K为0，数量占比为0.4，即包含2个振子。各振子相位正弦值随时间变化情况如下图所示。观察到3个第一类合作振子实现同步，而2个第二类孤立振子自始至终按自身相位做周期运动。
其次测试孤立振子占优的情形。假设第一类合作振子数量占比p为0.4，即包含2个振子，第二类孤立振子数量占比为0.6，即包含3个振子。各振子相位正弦值随时间变化情况如下图所示。观察到3个第二类孤立振子自始至终按自身相位做周期运动，而2个第一类合作振子和孤立振子4同步。换言之，孤立振子的表现始终“我自岿然不动”，而合作振子此时“选择和某个孤立振子合作”。
总的来看，当系统由合作振子和孤立振子组成时，孤立振子不受平均场影响，始终按自身相位运动，合作振子之间形成同步。合作振子可能以单独的相位运动，也可能和某个孤立振子同步运动，具体为何种方式取决于振子的初始相位，这里不作展开讨论。
合作振子和延迟振子
下面我们测试合作振子和延迟振子组成的系统。延迟振子的耦合强度K为正数，受平均场吸引，但是和平均场存在β的相位延迟。
首先测试合作振子占优的情形。假设第一类合作振子的数量占比p为0.6，即包含3个振子，不存在相位延迟，β1＝0；第二类延迟振子耦合强度K为0.5，数量占比为0.4，即包含2个振子，相位延迟β2＝π/2。
各振子相位正弦值随时间变化情况如下图所示。观察到3个第一类合作振子聚集成团簇，从而实现同步，而2个第二类延迟振子聚集成另一个团簇，也实现同步。两个团簇的相位差为相位延迟β2＝π/2。
其次测试延迟振子占优的情形。假设第一类合作振子的数量占比p为0.4，即包含2个振子，不存在相位延迟，β1＝0；第二类延迟振子耦合强度K为0.5，数量占比为0.6，即包含3个振子，相位延迟β2＝π/2。
各振子相位正弦值随时间变化如下图所示。观察到系统最终状态与合作振子占优的情形一致，即合作振子和延迟振子分别形成团簇而实现同步，团簇的相位差为相位延迟β2。和上一种情形略有差异之处在于，由于合作振子数量更少，同步的速度相对更慢。
合作振子占比p的影响
通过上述测试，我们发现当系统由合作振子和叛逆振子组成时，系统最终状态和合作振子占比p密切相关。本节我们将进一步讨论合作振子占比p的影响。本节模拟振子数量N从5扩展到100，第二类叛逆振子的耦合强度K为-0.5，全部振子的自然频率取自半高宽γ为0.05的洛伦兹分布。
我们对合作振子占比p进行遍历，计算每种情况系统最后2000步模拟序参量的均值。如下图所示，序参量r均值随着p的增大而提升，即合作振子占比越高，系统越倾向于同步。
从上图还能观察到一个现象：序参量均值并非平滑地提升，而是在p＝0.4~0.6之间存在一处“突起”。这一现象在Hong－Strogatz模型的原始文献Kuramoto Model of Coupled Oscillators with Positive and Negative Coupling Parameters中也有体现。该文献中，Hong和Strogatz（2011）模拟了由25600个振子构成的系统，包含合作和叛逆两类振子，其它参数和我们模拟所使用的参数一致。该研究得到的不同合作振子占比下系统序参量均值如下图所示，同样能够观察到p＝0.4~0.6之间的“突起”。
“突起”背后有何含义？实际上代表了一种特殊的系统状态——行波态。在合作振子和叛逆振子组成的系统中，根据不同的合作振子占比p值，存在三种不同的系统状态：非相干态（Incoherent State）、π态（π-state）和行波态（Traveling Wave State）。下面我们分别举例说明。
非相干态
当合作振子占比p为0.1时，系统序参量随时间变化情况如左下图，序参量r接近0，即接近完全不同步的状态。之所以称“接近”是因为r不完全等于0。
我们分别取两类振子的最终相位，绘制其频次分布图，如右上图所示。两类振子的最终相位均匀分布在[0, 2π]的区间内，说明振子没有出现同步现象，即使是合作振子也未同步。这种状态称为非相干态。
π态
当合作振子占比p为0.7时，系统序参量随时间变化情况如左下图，序参量r在0.4附近波动，表现出部分同步。两类振子最终相位频次分布如右下图所示。两类振子各自聚集成团簇，如同此前提到的“物以类聚，人以群分”。
此时，合作振子相位均值为4.81，叛逆振子相位中位数为1.88，两类振子相位差为2.9，接近π≈3.14。考虑到自然频率洛伦兹分布采样的随机性，以及数值模拟带来的误差，该相位差理论上应等于π，即两个团簇相位恰好相反。这种状态称为π态，是一种相对稳定的状态。
行波态
当合作振子占比p为0.55时，系统序参量随时间变化情况如左下图，序参量r在0.1~0.4的范围内波动，表现出部分同步。两类振子最终相位频次分布如右下图所示。两类振子各自聚集成团簇，其中合作振子的团簇更紧密，分布更窄，而叛逆振子的团簇较松散，分布更宽。
此时，合作振子相位均值为2.34，叛逆振子相位中位数为4.08，两类振子相位差绝对值为1.7，接近π/2的水平，离π较远。这种状态称为行波态。行波态和π态的区别在于，在旋转坐标系下，π态系统平均相位不随时间变化，而行波态系统平均相位以固定的速度运动，故称之为“行波”。关于藏本模型π态和行波态的判断涉及更为复杂的理论推导，本文不作更深入的展开。
小结
本章我们对Sakaguchi－Kuramoto模型进行模拟，探讨振子合作、叛逆、孤立、延迟等行为以及合作振子数量占比，对系统最终状态的影响，核心结论如下：
1. 当系统由合作振子和叛逆振子组成时：若合作振子占优，则两类振子分别形成团簇而实现同步，团簇间相位差等于π，即处于反相位状态；若叛逆振子占优，则系统无法实现同步。
2. 当系统由合作振子和孤立振子组成时，孤立振子不受平均场影响，始终按自身相位运动，合作振子间形成同步。
3. 当系统由合作振子和延迟振子组成时，两类振子分别形成团簇而实现同步，团簇间相位差等于相位延迟。
4. 合作振子占比越高，系统越倾向于同步；根据不同的合作振子占比p值，存在非相干态、π态和行波态三种不同的系统状态。
总结
人们常说生命是奇迹。动物心脏的大量起搏细胞以同一种节律放电，合奏出生命的不朽乐章。鸟鸣蝉噪，鱼贯雁行，流萤闪烁，落英缤纷，大自然并不存在一位全知全能的指挥者，而万物却自发地以相同的频率鸣唱，以相同的步调行进，以相同的节奏经历着成熟和衰亡。人们也常说文明是奇迹。中国工人和北美农民远隔重洋素不相识，却共同用劳动书写经济的高度繁荣，也在危机来临之际环球共此冷暖。
当我们探索奇迹背后的奥秘，我们发现奇迹并不是造物主提前设计好的剧本，而是参与其中的每一个细胞、个体、企业自发地同步所产生的合力。系统的每个微观个体都具备一定周期性，众多个体通过相互间的耦合，最终形成步调一致的宏观周期。这种现象在物理学中称为耦合振子同步。而藏本模型是刻画耦合振子同步现象的经典模型。
藏本模型的基本形式是描述系统中每个振子相位的一维动力学方程。模型中振子相位的变化量（即运动速度）由1）该振子本身自然频率ω，2）耦合强度K，3）该振子和其它振子的相位差决定。耦合强度K＞0使得振子间相互吸引。引入系统序参量r和平均相位ψ和平均场的概念后，振子受到其它振子的吸引而趋于同步，等价于振子受到平均场的吸引而趋于同步。
藏本模型中振子数量N、耦合强度K、自然频率洛伦兹分布半高宽γ这三项关键参数对同步的影响不尽相同。当振子数量N≥3时，振子数量对同步速度影响不大。耦合强度K小于某个定值时，无法同步；大于某个定值时，K越大同步越快。振子自然频率服从洛伦兹分布的半高宽γ越大，系统越难同步。当振子的自然频率介于-Kr和Kr之间时，振子能与平均场同步；当振子的自然频率超出该范围时，振子无法同步。
Sakaguchi－Kuramoto模型是藏本模型的变式，可借助该模型研究振子合作、叛逆等行为对系统终态的影响。当系统由合作和叛逆振子组成时：若合作振子占优，则两类振子分别形成团簇而同步，团簇间相位差为π，即处于反相位；若叛逆振子占优，则系统无法同步。当系统由合作和孤立振子组成时，孤立振子始终按自身相位运动，合作振子形成同步。当系统由合作和延迟振子组成时，两类振子分别形成团簇而同步，团簇间存在相位延迟。合作振子占比越高，系统越倾向于同步；根据不同的合作振子占比p值，存在非相干态、π态和行波态三种不同的系统状态。
藏本模型目前主要运用在自然科学领域。事实上，在人类的经济活动和金融系统中也存在类似的耦合振子同步现象。例如多数企业的库存具有周期性，这是企业应对外部不确定性时的缓冲机制；商品、货币和信息的流通使得企业与企业之间相互影响和制约；微观企业相互耦合逐渐同步，最终形成宏观级别的周期。藏本模型如何应用于经济周期的刻画，是否可能有助于我们理解金融经济系统周期的起源，这些问题值得进一步探索。
参考文献
鞠萍. (2016). 耦合振子系统中的动力学行为. (Doctoral dissertation).
Hong, H. , & Strogatz, S. H. . (2011). Kuramoto model of coupled oscillators with positive and negative coupling parameters: an example of conformist and contrarian oscillators. Physical Review Letters, 106(5), p.054102.1-054102.4.
Neda, Z. , Ravasz, E. , Brechet, Y. , Vicsek, T. , & Barabasi, A. L. . (2000). Physics of the rhythmic applause. Physical Review E Statistical Physics Plasmas Fluids & Related Interdisciplinary Topics, 61(6).
Neda, Z. , Ravasz, E. , Brechet, Y. , Vicsek, T. , & Barabasi, A. L. . (2000). The sound of many hands clapping - tumultuous applause can transform itself into waves of synchronized clapping. Nature, 403(6772), 849-850.
风险提示
包括藏本模型在内的动力学系统模型是对真实世界的简单刻画，存在过度简化的可能；周期规律基于历史数据总结，历史规律可能失效；周期规律对市场长期规律进行研究，无法判断短期的市场情绪与政策冲击带来的影响；市场在某些极端情形下可能出现规律以外的交易风险。
附录
藏本模型两种形式恒等变换的推导
本节展示藏本模型两种形式恒等变换的具体推导过程。
藏本模型静止坐标系转换至旋转坐标系的推导
本节展示藏本模型静止坐标系转换至旋转坐标系的具体推导过程。
免责声明与评级说明
公众平台免责申明
本公众平台不是华泰证券研究所官方订阅平台。相关观点或信息请以华泰证券官方公众平台为准。根据《证券期货投资者适当性管理办法》的相关要求，本公众号内容仅面向华泰证券客户中的专业投资者，请勿对本公众号内容进行任何形式的转发。若您并非华泰证券客户中的专业投资者，请取消关注本公众号，不再订阅、接收或使用本公众号中的内容。因本公众号难以设置访问权限，若给您造成不便，烦请谅解！本公众号旨在沟通研究信息，交流研究经验，华泰证券不因任何订阅本公众号的行为而将订阅者视为华泰证券的客户。
本公众号研究报告有关内容摘编自已经发布的研究报告的，若因对报告的摘编而产生歧义，应以报告发布当日的完整内容为准。如需了解详细内容，请具体参见华泰证券所发布的完整版报告。
本公众号内容基于作者认为可靠的、已公开的信息编制，但作者对该等信息的准确性及完整性不作任何保证，也不对证券价格的涨跌或市场走势作确定性判断。本公众号所载的意见、评估及预测仅反映发布当日的观点和判断。在不同时期，华泰证券可能会发出与本公众号所载意见、评估及预测不一致的研究报告。
在任何情况下，本公众号中的信息或所表述的意见均不构成对客户私人投资建议。订阅人不应单独依靠本订阅号中的信息而取代自身独立的判断，应自主做出投资决策并自行承担投资风险。普通投资者若使用本资料，有可能会因缺乏解读服务而对内容产生理解上的歧义，进而造成投资损失。对依据或者使用本公众号内容所造成的一切后果，华泰证券及作者均不承担任何法律责任。
本公众号版权仅为华泰证券股份有限公司所有，未经公司书面许可，任何机构或个人不得以翻版、复制、发表、引用或再次分发他人等任何形式侵犯本公众号发布的所有内容的版权。如因侵权行为给华泰证券造成任何直接或间接的损失，华泰证券保留追究一切法律责任的权利。本公司具有中国证监会核准的“证券投资咨询”业务资格，经营许可证编号为：91320000704041011J。
华泰金工深度报告一览
金融周期系列研究（资产配置）
【华泰金工林晓明团队】2020年中国市场量化资产配置年度观点——周期归来、机会重生，顾短也兼长20200121
【华泰金工林晓明团队】量化资产配置2020年度观点——小周期争明日，大周期赢未来20200116
【华泰金工林晓明团队】风险预算模型如何度量风险更有效-改进风险度量方式稳定提升风险模型表现的方法
【华泰金工林晓明团队】周期双底存不确定性宜防守待趋势——短周期底部拐头机会渐增，待趋势明朗把握或更大20191022
【华泰金工林晓明团队】二十年一轮回的黄金投资大周期——黄金的三周期定价逻辑与组合配置、投资机会分析20190826
【华泰金工林晓明团队】如何有效判断真正的周期拐点？——定量测度实际周期长度提升市场拐点判准概率
【华泰金工林晓明团队】基钦周期的长度会缩短吗？——20190506
【华泰金工林晓明团队】二十载昔日重现，三四年周期轮回——2019年中国与全球市场量化资产配置年度观点（下）
【华泰金工林晓明团队】二十载昔日重现，三四年周期轮回——2019年中国与全球市场量化资产配置年度观点（上）
【华泰金工林晓明团队】周期轮动下的BL资产配置策略
【华泰金工林晓明团队】周期理论与机器学习资产收益预测——华泰金工市场周期与资产配置研究
【华泰金工林晓明团队】市场拐点的判断方法
【华泰金工林晓明团队】2018中国与全球市场的机会、风险 · 年度策略报告（上）
【华泰金工林晓明团队】基钦周期的量化测度与历史规律 · 华泰金工周期系列研究
【华泰金工林晓明团队】周期三因子定价与资产配置模型（四）——华泰金工周期系列研究
【华泰金工林晓明团队】周期三因子定价与资产配置模型（三）——华泰金工周期系列研究
【华泰金工林晓明团队】周期三因子定价与资产配置模型（二）——华泰金工周期系列研究
【华泰金工林晓明团队】周期三因子定价与资产配置模型（一）——华泰金工周期系列研究
【华泰金工林晓明团队】华泰金工周期研究系列 · 基于DDM模型的板块轮动探索
【华泰金工林晓明团队】市场周期的量化分解
【华泰金工林晓明团队】周期研究对大类资产的预测观点
【华泰金工林晓明团队】金融经济系统周期的确定（下）——华泰金工周期系列研究
【华泰金工林晓明团队】金融经济系统周期的确定（上）——华泰金工周期系列研究
【华泰金工林晓明团队】全球多市场择时配置初探——华泰周期择时研究系列
行业指数频谱分析及配置模型：市场的周期分析系列之三
【华泰金工林晓明团队】市场的频率——市场轮回，周期重生
【华泰金工林晓明团队】市场的轮回——金融市场周期与经济周期关系初探
周期起源
【华泰金工林晓明团队】周期在供应链管理模型的实证——华泰周期起源系列研究之六
【华泰金工林晓明团队】不确定性与缓冲机制——华泰周期起源系列研究报告之五
【华泰金工林晓明团队】周期是矛盾双方稳定共存的结果——华泰周期起源系列研究之四
【华泰金工林晓明团队】周期是不确定性条件下的稳态——华泰周期起源系列研究之三
【华泰金工林晓明团队】周期趋同现象的动力学系统模型——华泰周期起源系列研究之二
【华泰金工林晓明团队】从微观同步到宏观周期——华泰周期起源系列研究报告之一
