Adoption of deep learning Markov model combined with copula function in portfolio risk measurement

The model based on Copula function can be adopted to study the correlation between variables and their characteristics. The binary Copula function is defined as follows.

The common Archimedes Copula functions include the Gumbel Copula function, the Clayton Copula function and the Frank Copula function [14, 15]. The mathematical expressions of these three functions are as follows, respectively. (3) C(u1,u2,⋯,ud;θ1)=exp{−[∑k=1d(−lnuk)1θ1]θ1} C\left( {{u_1},{u_2}, \cdots ,{u_d};{\theta _1}} \right) = \exp \left\{ { - {{\left[ {\sum\limits_{k = 1}^d {{\left( { - \ln {u_k}} \right)}^{{1 \over {{\theta _1}}}}}} \right]}^{{\theta _1}}}} \right\} (4) C(u1,u2,⋯,ud;θ2)=(∑k=1duk−θ2−d+1)−1θ2 C\left( {{u_1},{u_2}, \cdots ,{u_d};{\theta _2}} \right) = {\left( {\sum\limits_{k = 1}^d u_{_k}^{ - {\theta _2}} - d + 1} \right)^{ - {1 \over {{\theta _2}}}}} (5) C(u1,u2,⋯,ud;θ3)=1θ3ln(1+∏k=1d(e−θ3uk−1)(e−θ3−1)d−1) C\left( {{u_1},{u_2}, \cdots ,{u_d};{\theta _3}} \right) = {1 \over {{\theta _3}}}\ln \left( {1 + {{\prod\limits_{k = 1}^d \left( {{e^{ - {\theta _3}{u_k}}} - 1} \right)} \over {{{\left( {{e^{ - {\theta _3}}} - 1} \right)}^{d - 1}}}}} \right) In Eq. (3), (−lnu)1θ1 {\left( { - \ln u} \right)^{{1 \over {{\theta _1}}}}} is the generator of Gumbel Copula function, θ₁ ∈ (0, 1], and the closer it approaches to 0, the closer it approaches to complete correlation; in Eq. (4), u^−θ₂ − 1 is the generator of Clayton Copula function, θ₂ ∈ (0, ∞), and the closer it approaches to ∞, the closer it approaches to complete correlation; lne−θ3u−1e−θ3−1 \ln {{{e^{ - {\theta _3}u}} - 1} \over {{e^{ - {\theta _3}}} - 1}} is the generator of the Frank Copula function, θ₃ ∈ (0, ∞), and the closer it approaches to ∞, the closer it approaches to complete correlation.

Then the Gumbel Copula, Clayton Copula and Frank Copula functions are combined to build a mixed Copula function, and the weight values of these three functions are w₁, w₂, w₃. The mixed Copula function can be expressed as follows. (6) Cm(u1,u2,⋯,ud;ψ)=∑k=13wkCk(u1,u2,⋯,ud;θk) {C_m}\left( {{u_1},{u_2}, \cdots ,{u_d};\psi } \right) = \sum\limits_{k = 1}^3 {w_k}{C_k}\left( {{u_1},{u_2}, \cdots ,{u_d};{\theta _k}} \right) In the above equation, 0 ≤ w_k ≤ 1, and are all uniform variables in [0,1]; ψ = {θ₁, θ₂, θ₃, w₁, w₂, w₃} is the parameter set.

Then the density function of the mixed Copula function is as follows. (7) Cm(u1,u2,⋯,ud;ψ)=∂dCm(u1,u2,⋯,ud;ψ)∂u1∂u2⋯∂ud=∑k=13wkCk(u1,u2,⋯,ud;θk) {C_m}\left( {{u_1},{u_2}, \cdots ,{u_d};\psi } \right) = {{{\partial ^d}{C_m}\left( {{u_1},{u_2}, \cdots ,{u_d};\psi } \right)} \over {\partial {u_1}\partial {u_2} \cdots \partial {u_d}}} = \sum\limits_{k = 1}^3 {w_k}{C_k}\left( {{u_1},{u_2}, \cdots ,{u_d};{\theta _k}} \right)

The state in HMM is hidden, which can be obtained by the observation sequence and transition probability matrix between states. The four parameters in HMM can be described as the number of states, distribution of initial probability, transition probability matrix and density function of observation sequence [16]. It is assumed that there are HM series {Xh} with m states, and d dimensional continuous observation series {Yh} has the conditional probability density function P_i (y_h) = P(Y_h = y_h |X_h = i) when the state is i. Therefore, the unconditional probability vector of the initial probability distribution is as follows. (8) π1*m=[P(X1=1),P(X1=2),⋯,P(X1=m)] {\pi _{1*m}} = \left[ {P\left( {{X_1} = 1} \right),P\left( {{X_1} = 2} \right), \cdots ,P\left( {{X_1} = m} \right)} \right] Assuming that the observed sequence is {Yh} and the state sequence is {Xh}, the hidden Markov chain needs to meet the following conditions.

According to the above two properties, the theory of conditional probability is combined with the theory of probability in mathematical statistics to carry out the derivation of P(Y1,⋯,Yh,X1,⋯,Xh): P(Y1,⋯,Yh,X1,⋯,Xh)=∏h=1Hpxh(yh)∏h=1H−1axxh+1πx1. P\left( {{Y_1}, \cdots ,{Y_h},{X_1}, \cdots ,{X_h}} \right):\;P\left( {{Y_1}, \cdots ,{Y_h},{X_1}, \cdots ,{X_h}} \right) = \prod\limits_{h = 1}^H {p_{{x_h}}}\left( {{y_h}} \right)\prod\limits_{h = 1}^{H - 1} {a_{x{x_{h + 1}}}}{\pi _{{x_1}}}. Based on the HMM model, the Copula model is nested in the framework to obtain the dynamic hybrid Copula model. The construction process is as follows.

Assuming that the matrix of order H × D is Y, the matrix formed by the d dimensional observation sequence at time h is {Yh}; the matrix of order H × m order is X, the matrix formed by m state sequences at time h is {Xh}; the matrix of order H × n order is K, the matrix formed by the n-dimensional mixed Copula function branch sequence at the time h is {Rh}.

According to Sklar theorem [17], the mathematical expression of the joint density function f_i of d dimension in the i state at time h is as follows. (9) fi(yh1,⋯,yhd)=ci(F1(yh1;η1)⋯,Fd(yhd;ηd);ψ)∏l=1dfl(yhl;ηl) {f_i}\left( {{y_{h1}}, \cdots ,{y_{hd}}} \right) = {c_i}\left( {{F_1}\left( {{y_{h1}};{\eta _1}} \right) \cdots ,{F_d}\left( {{y_{hd}};{\eta _d}} \right);\psi } \right)\prod\limits_{l = 1}^d {f_l}\left( {{y_{hl}};{\eta _l}} \right) In Eq. (9), c_i () is the mixed Copula density function of i state; η is the parameter vector of edge distribution function and ψ = {θ₁, θ₂, θ₃, w₁, w₂, w₃} is the parameter set.

If it is the parameter set, ζ = {π_i, a_ij, w_in, θ_in, η_l} is the dependent parameter of the nth function of i state, then the mathematical expression of the HMM likelihood function of d dimension is as follows. (10) L(ζ,Y,X,R)=P(Y,X,R|ζ)=πx1∏h=1H−1axhxh+1∏h=1H(∏l=1df(yhl;ηl)wxhKhcxhKh(uh)) L\left( {\zeta ,Y,X,R} \right) = P\left( {Y,X,R\left| \zeta \right.} \right) = {\pi _{{x_1}}}\prod\limits_{h = 1}^{H - 1} {a_{{x_h}{x_{h + 1}}}}\prod\limits_{h = 1}^H \left( {\prod\limits_{l = 1}^d f\left( {{y_{hl}};{\eta _l}} \right){w_{{x_h}{K_h}}}{c_{{x_h}{K_h}}}\left( {{u_h}} \right)} \right) Observation sequence and parameter value vector ζ are known, and then the mathematical expressions for the conditional expectation Q of the completely likelihood function that is at K step of iteration is as follows. (11) Q(ζ,ζ(r))=∑{X,R}P(Y,X,R|Y;ζ(r))logP(Y,X,R|ζ)=∑i=1mlogπiγi(k)(1)+∑i=1m∑j=1m∑h=1H−1logaijεij(k)(h)+ ∑i=1m∑h=1Hγi(k)(h)log(∏l=1dfl(yhl;ηl))+ ∑i=1m∑h=1H∑v=13γi(k)(t)(logciv(uh,θiv|ψ(k))+log(wiv|ψ(k)))wiv∑v=13wivciv \matrix{ {Q\left( {\zeta ,{\zeta _{\left( r \right)}}} \right)} \hfill & { = \sum\limits_{\left\{ {X,R} \right\}} P\left( {Y,X,R\left| {Y;{\zeta _{\left( r \right)}}} \right.} \right)\log P\left( {Y,X,R\left| \zeta \right.} \right) = \sum\limits_{i = 1}^m \log {\pi _i}\gamma _i^{\left( k \right)}\left( 1 \right) + \sum\limits_{i = 1}^m \sum\limits_{j = 1}^m \sum\limits_{h = 1}^{H - 1} \log {a_{ij}}\varepsilon _{ij}^{\left( k \right)}\left( h \right)} \hfill \cr {} \hfill & { + \;\sum\limits_{i = 1}^m \sum\limits_{h = 1}^H \gamma _i^{\left( k \right)}\left( h \right)\log \left( {\prod\limits_{l = 1}^d {f_l}\left( {{y_{hl}};{\eta _l}} \right)} \right)} \hfill \cr {} \hfill & { + \;\sum\limits_{i = 1}^m \sum\limits_{h = 1}^H \sum\limits_{v = 1}^3 \gamma _i^{\left( k \right)}\left( t \right)\left( {\log {c_{iv}}\left( {{u_h},{\theta _{iv}}\left| {{\psi ^{\left( k \right)}}} \right.} \right) + \log \left( {{w_{iv}}\left| {{\psi ^{\left( k \right)}}} \right.} \right)} \right){{{w_{iv}}} \over {\sum\limits_{v = 1}^3 {w_{iv}}{c_{iv}}}}} \hfill \cr } In the above equation, the conditional expectation Q function consists of five parts. Therefore, it needs to iterate to step K to maximise the Q function, which needs K steps to maximise each part. If the weights of the initial distribution, transition matrix and Copula functions are conditional, the Lagrange multiplication can be adopted to maximise the parameters and then the maximum likelihood function is adopted to estimate the edge distribution. In Eq. (11), the dependent parameters are multivariate nonlinear mathematical expressions with boundary conditions. Taking the optimisation of Autoencoder model in deep learning as an example [18], the L-BFGS can solve similar problems with gradient descent function, SGD and other functions, as shown in Figure 1. But in most cases, the L-BFGS function converges faster and has less memory overhead than other algorithms [19].

Fig. 1

The L-BFGS-B algorithm is adopted to obtain the estimated value after k+1 iteration, and the mathematical expression of the final estimated value is as follows. (12) πi=γi(k)(1) {\pi _i} = \gamma _i^{\left( k \right)}\left( 1 \right) (13) aij=∑h=1H−1εij(k)(h)∑j=1m∑h=1H−1εij(k)(h) {a_{ij}} = {{\sum\limits_{h = 1}^{H - 1} \varepsilon _{ij}^{\left( k \right)}\left( h \right)} \over {\sum\limits_{j = 1}^m \sum\limits_{h = 1}^{H - 1} \varepsilon _{ij}^{\left( k \right)}\left( h \right)}} (14) wiv=∑h=1Hγi(k)(t)wiv(k)civ(uh,θiv(k))∑v=13wiv(k)civ(uh,θiv(k))∑h=1Hγi(k)(t) {w_{iv}} = {{\sum\limits_{h = 1}^H \gamma _i^{\left( k \right)}\left( t \right){{w_{iv}^{\left( k \right)}{c_{iv}}\left( {{u_h},\theta _{iv}^{\left( k \right)}} \right)} \over {\sum\limits_{v = 1}^3 w_{iv}^{\left( k \right)}{c_{iv}}\left( {{u_h},\theta _{iv}^{\left( k \right)}} \right)}}} \over {\sum\limits_{h = 1}^H \gamma _i^{\left( k \right)}\left( t \right)}}

To construct the portfolio model with Copula function, it is necessary to determine the edge distribution of Copula and determine the edge distribution function. The appropriate Copula function is selected to describe the dependence structure relation of multi-variable edge distribution. The generalised autoregressive conditional heteroskedast (GARCH) model is adopted to generalise the marginal distribution of variables [20]. Copula function is adopted to combine edge distribution functions of different variables, and the resulting joint distribution function is the portfolio function of financial assets. Assuming that the financial time series are {R1t} and {R2t}, the mathematical expression of GARCH (1,1) model of the financial time series is as follows. (15) {rit=μi+εitεit=ηitσit, ηit∼i,t,d Fi()σit2=ωi+αiεit−12+βiσit−12(η1t,η2t)|Ft−1∼Ct(F1(η1t),F2(η2t)|Ft−1) \left\{ {\matrix{ {{r_{it}} = {\mu _i} + {\varepsilon _{it}}} \hfill \cr {{\varepsilon _{it}} = {\eta _{it}}{\sigma _{it}},\quad {\eta _{it}}\mathop \sim \limits^{{\rm i,t,d}{\kern 1pt} } {F_i}\left( {} \right)} \hfill \cr {\sigma _{_{it}}^2 = {\omega _i} + {\alpha _i}\varepsilon _{_{it - 1}}^2 + {\beta _i}\sigma _{_{it - 1}}^2} \hfill \cr {\left( {{\eta _{1t}},{\eta _{2t}}} \right)\left| {{F_{t - 1}} \sim {C_t}\left( {{F_1}\left( {{\eta _{1t}}} \right),{F_2}\left( {{\eta _{2t}}} \right)\left| {{F_{t - 1}}} \right.} \right)} \right.} \hfill \cr } } \right. where r_it is the yield rate, η_it obeys F_i() distribution, α is the coefficient of ARCH and β is the coefficient of GARCH.

It contains binary normal Copula and T-Copula models, whose mathematical expressions are as follows. (16) C(u,v;p)=∫−∞φ−1(u)∫−∞φ−1(v)12π1−p2exp[−s2−2pst+t22(1−p2)]dsdt C\left( {u,v;p} \right) = \int_{ - \infty }^{{\varphi ^{ - 1}}\left( u \right)} \int_{ - \infty }^{{\varphi ^{ - 1}}\left( v \right)} {1 \over {2\pi \sqrt {1 - {p^2}} }}\exp \left[ { - {{{s^2} - 2pst + {t^2}} \over {2\left( {1 - {p^2}} \right)}}} \right]dsdt (17) C(u,v;p,k)=∫−∞tk−1(u)∫−∞tk−1(v)12π1−p2[1+s2−2pst+t2k(1−p2)]−(k+2)2dsdt C\left( {u,v;p,k} \right) = \int_{ - \infty }^{t_k^{ - 1}\left( u \right)} \int_{ - \infty }^{t_k^{ - 1}\left( v \right)} {1 \over {2\pi \sqrt {1 - {p^2}} }}{\left[ {1 + {{{s^2} - 2pst + {t^2}} \over {k\left( {1 - {p^2}} \right)}}} \right]^{{{ - \left( {k + 2} \right)} \over 2}}}dsdt Then Monte Carlo simulation is adopted to simulate the value at risk (VaR) of the investment portfolio [21]. Assuming that the initial total investment is fixed and the time range is demarcated as 1 day, the total assets of the X and Y portfolios of the two types of investments at time t is V_t = n₁P_1,t + n₂P_2,t, where n is the proportion of these two types of investments in the total portfolio, and P is the value of these two types of investments at time t. And the logarithmic yield on the portfolio is as follows. (18) Rt+1=log(n1P1,tn1P1,t+n2P2,teX+n1P1,tn1P1,t+n2P2,teY)=log(λ1eX+λ2eY) {R_{t + 1}} = \log \left( {{{{n_1}{P_{1,t}}} \over {{n_1}{P_{1,t}} + {n_2}{P_{2,t}}}}{e^X} + {{{n_1}{P_{1,t}}} \over {{n_1}{P_{1,t}} + {n_2}{P_{2,t}}}}{e^Y}} \right) = \log \left( {{\lambda _1}{e^X} + {\lambda _2}{e^Y}} \right) λ is the weight of the two types of investment in the portfolio, the sum of which is 1.

VaR satisfies the following equation. (19) 1−α=P(R≤−r)=∫−∞−r−log(−λt)dx∫−∞log(e−rλ2−λ1λ2ey)ctft(x|Ft−1)gt(x|Ft−1)dy 1 - \alpha = P\left( {R \le - r} \right) = \int_{ - \infty }^{ - r - \log \left( { - {\lambda _t}} \right)} dx\int_{ - \infty }^{\log \left( {{e^{ - r}}{\lambda _2} - {\lambda _1}{\lambda _2}{e^y}} \right)} {c_t}{f_t}\left( {x\left| {{F_{t - 1}}} \right.} \right){g_t}\left( {x\left| {{F_{t - 1}}} \right.} \right)dy α is the confidence level.

In practice, the solution of Eq. (19) is more difficult, so Monte Carlo simulation can be considered for solution simulation; the simulation steps are shown in Figure 2.

Fig. 2

Recurrent neural network (RNN) is used for prediction; RNN is also known as a deep learning Markov model, as its properties conform to the specific generalised Markov model [22, 23]. The basic structure of RNN model is shown in Figure 3.

Fig. 3

The neurons in the hidden layer in the RNN are affected by the neurons in the upper layer. When t = 0, there is no state of the hidden layer neurons at the previous moment. Different from other models, the parameters in the RNN model are shared with each other, which makes the RNN model carry out the same processing at different moments, which reduces the workload of parameter learning for the model [24]. At time t, the gradient of the error function E in the cyclic layer weight matrix W with respect to RNN is as follows. (20) ∂E∂Wt=∂E∂rt∂rt∂Wt=[δ1th1t−1δ1th2t−1⋯δ1thnt−1δ2th1t−1δ2th2t−1⋯δ2thnt−1⋮⋮⋮⋮δnth1t−1δnth2t−1⋯δnthnt−1] {{\partial E} \over {\partial {W_t}}} = {{\partial E} \over {\partial {r_t}}}{{\partial {r_t}} \over {\partial {W_t}}} = \left[ {\matrix{ {\delta _1^th_1^{t - 1}} & {\delta _1^th_2^{t - 1}} & \cdots & {\delta _1^th_n^{t - 1}} \cr {\delta _2^th_1^{t - 1}} & {\delta _2^th_2^{t - 1}} & \cdots & {\delta _2^th_n^{t - 1}} \cr \vdots & \vdots & \vdots & \vdots \cr {\delta _n^th_1^{t - 1}} & {\delta _n^th_2^{t - 1}} & \cdots & {\delta _n^th_n^{t - 1}} \cr } } \right] Then the gradient of the error function E in the weight matrix W in the cyclic layer is the sum ∂E∂W=∑i=1t∂E∂Wt {{\partial E} \over {\partial W}} = \sum\limits_{i = 1}^t {{\partial E} \over {\partial {W_t}}} of the gradients at each moment. Similarly, the gradient of the error function E to the weight matrix U can be obtained as ∂E∂U=∑i=1t∂E∂Ut {{\partial E} \over {\partial U}} = \sum\limits_{i = 1}^t {{\partial E} \over {\partial {U_t}}} .

After the parameter gradient of weight matrix in RNN is determined, the stochastic gradient descent algorithm is adopted to train RNN. Then the RNN prediction results are evaluated by calculating the error rate and the average error rate. The calculation equations of error rate e and average error rate at time t are as follows. (21) e=|y′−y|y*100% e = {{\left| {{y^\prime} - y} \right|} \over y}*100\% (22) e¯=1N∑t=1N|y′−y|y*100% \overline e = {1 \over N}\sum\limits_{t = 1}^N {{\left| {{y^\prime} - y} \right|} \over y}*100\% Y is the prediction results, y is the actual results and N is the data amount.

The bank (A), insurance (B), securities (C) and trust (D) in the SWS index (http://www.swsindex.com/) are taken as research objects. The closing prices of solstice from 2 January 2015 to 31 December 2019 are collected as the research object. As all of the data are time series, data preprocessing is required to ensure the stability of data before building the model by these data [25]. In this study, the logarithmic rate of yield r of various industries is selected for empirical analysis. The descriptive statistics of data of various industries are conducted, and the results are shown in Figure 4. Compared with the standard normal distribution, the distribution of exponential logarithmic yields of various industries has a more obvious kurtosis and skewness.

Fig. 4

Through descriptive statistics, the kurtosis values of logarithmic yields of these four industries series index are, respectively, 9.668, 6.268, 6.842 and 6.557, which all exceed the peak value 3 of normal distribution. As they have high peak value, they do not conform to the normal distribution, and the yield rates have the characteristics of peak distribution. According to the analysis of skewness value, the skewness value of the exponential logarithmic yield of these four industries is, respectively, 0.031, 0.082, 0.043 and −0.741, which indicates that except for D, all sequences A, B and C are right-biased. The P values are all <0.05 after Jarque Bera statistical test, indicating that they do not conform to the initial hypothesis of normal distribution.

The GARCH model is adopted to construct the exponential logarithmic yield series models of bank (A), insurance (B), securities (C) and trust (D). The final constructed model is the edge distribution function. After multiple verifications, it is found that GARCH (1,1) is the most effective method to verify all the yield rate sequence models. Therefore, normal distribution and T distribution are used in this study to simulate the residuals of yield series. The GARCH (1,1) models of normal distribution and t-distribution are constructed by MATLAB. The parameters of each model are estimated respectively. The final estimation results of edge distribution are shown in Table 2. Akaike information criterion (AIC) value can reflect the optimal benign of statistical model fitting. The estimation results of edge distribution parameters obtained from Table 1 show that the AIC value of normal distribution is greater than that of T-distribution, which indicates the fitting effect of edge distribution of exponential logarithmic yield rates of A, B, C and D industry sequences based on T distribution is obviously better than that of normal distribution. As A, B, C, and D exponential logarithmic yield sequence have the feature of thick tail, it indicates that the T distribution has an obvious advantage in describing the yield sequence with peak and thick tail distribution feature.

In order to eliminate autocorrelation, autoregressive moving average (ARMA) model [26] is used for data processing in this study. Compared with the GARCH model, the generalised autoregressive score (GAS) model can make full use of the density function, so it is adopted to overcome the heteroscedasticity in the data [27]. Assuming the yield is y, and y = u + σε, where u is the mean value equation, which satisfies ARMA model; σ is the conditional standard deviation, and the fluctuation in f = σ² satisfies GAS. The density function of y is p(y|f, F), where F is the information set. Then the mathematical expression of GAS (1,1) model is as follows. (23) {ft+1=w+αst+βftst=St(ft)∇t(yt,ft)St=Et−1[∇t,∇t′]−γ∇t=∂lnp(yt|ft,Ft)∂ft \left\{ {\matrix{ {{f_{t + 1}} = w + \alpha {s_t} + \beta {f_t}} \hfill\cr {{s_t} = {S_t}\left( {{f_t}} \right){\nabla _t}\left( {{y_t},{f_t}} \right)} \hfill\cr {{S_t} = {E_{t - 1}}{{\left[ {{\nabla _t},\nabla _{_t}^\prime} \right]}^{ - \gamma }}} \hfill\cr {{\nabla _t} = {{\partial \ln p\left( {{y_t}\left| {{f_t},{F_t}} \right.} \right)} \over {\partial {f_t}}}} \hfill\cr } } \right. In Eq. (23), ∇_t is the conditional score of the density function; s_t is the scale score of the observation sequence; S_t is the scalar; w, α, β are constant coefficients, the recursive direction of f is controlled by α, β < 1 and the persistence of recursion is controlled by it as γ = 0.

The GAS model is adopted to determine the sequence edge distribution of each industry, and the parameter estimation value of the GAS model is finally obtained, as shown in Table 3. As shown in Table 3 and Figure 5A, 5B and 5C all obey the Laplace distribution, which confirms that the sequences of the three industries have certain skewness, and it is the same as the analysis results of the GARCH model. The distribution of D industry is different from that of other industries, and it belongs to partial t-distribution. As shown in Table 3, the α values of these four industry sequences are all positive, which means that these sequences recourse in the fastest ascending algorithm. However, β > 0.95, which is close to 1, indicating that the fluctuation process of these four industry sequences has a long duration.

Fig. 5

Furthermore, the histogram results of the probability integral transform (PIT) values of the edge distribution in the 90% confidence interval of these four industry sequences are calculated. As shown in Figure 6, in the transformation process, the P value of A (banking) is manifested as P < 0.05 at the 3rd and 4th moments, while the P value of C (securities industry) is manifested as P < 0.05 at the final histogram. And the final KS test results (P values) of PIT values of all industries are >0.05, which indicates that the edge distributions fitted in this study also obey uniform distribution after PIT.

Fig. 6

The uniformly distributed variables obtained from the above PIT are calculated in the mathematical expressions of three single Gumbel Copula, Clayton Copula, Frank Copula and mixed Copula functions, which are Eqs 3–6. The parameters of the Gumbel Copula, Clayton Copula, Frank Copula, and mixed Copula functions are estimated, respectively, and the iteration results of parameters in different functions are shown in Figure 7. The expectation-maximisation (EM) algorithm is adopted for parameter iteration, and the parameters of different functions gradually become stable after 10 times of iteration.

Fig. 7

Therefore, the parameter estimates of the three single Copula functions are taken as the initial values of the parameters of the mixed Copula function, and the obtained parameter estimates of the single Copula and mixed Copula functions are shown in Table 4. According to the principle of maximum log-likelihood function and minimum AIC and Bayesian information criterion (BIC) [28], the log-likelihood function of the mixed Copula function is the largest and its AIC and BIC values are the smallest, which indicates the mixed Copula model constructed in this study has better effects compared with the single Copula model.

As shown in Figure 8, the relevant parameters in the HMM-based mixed Copula model gradually become stable after 50 iterations, and the parameters of the HMM-based mixed Copula model in the final state 1 are θ₁ = 1.208, θ₂ = 0.923, θ₃ = 8.926, w₁ = 0.505, w₂ = 0.482, w₃ = 0.025 and initial distribution π = 0.026. The parameters of the HMM-based mixed Copula model in the final state 2 are θ₁ = 2.153, θ₂ = 2.101, θ₃ = 18.220, w₁ = 0.553, w₂ = 0.339, w₃ = 0.087 and initial distribution π = 0.977.

Fig. 8

As shown in Table 4, the log-likelihood function of the mixed Copula model based on HMM is 2059.63, AIC = −4088.39, and BIC = −4010.98. State 1 is a low dependent state, while state 2 is a high dependent state. After analysing the probability of Copula model based on HMM in these two states, the probability of state 2 is the highest (97.55%), which indicates that these four industries are in a highly dependent state. As shown in Table 4, the mixed Copula model based on HMM proposed in this study has the largest logarithmic likelihood function and the lowest AIC and BIC values, which indicates that the model effect of nesting mixed Copula functions into HMM is optimal, which is consistent with the research results of Chun et al. (2015) [29].

Subsequently, the dynamic transition diagrams of states 1 and 2 of the sample from 2015 to 2019 are obtained, as shown in Figure 9A. These two dependent states will change at a certain frequency, yet the evaluation of this transformation is not regular and mainly depends on the degree of interdependence between industries. Studies show that the increase in interdependence will lead to the occurrence of risk infection, and the probability of risk infection is very high during this period. As shown in the low dependence probability graph (Figure 9B) and the high dependence probability graph (Figure 9C), the probability graph is basically consistent with the trend of the dynamic transition graph of states. Therefore, when an industry is in state 2 (high dependence state), there is a very strong dependence relationship between industries, which may cause macro or systemic risks from the perspective of dependency and risk contagion [30].

Fig. 9

The interdependence of industries under these two states is calculated, and the results are shown in Table 5. Under different states, the tail dependence between A (banking) and B (insurance) is the strongest, which indicates that they have the most sensitivity to the impact of risks, and are more susceptible to be impacted by the own risks of extreme events.

Assuming that the investor invests according to the weight of [0.5, 0.5] and the investment amount is 1 RMB. Monte Carlo simulation is adopted to measure the VaR value of various fields for 5,000 times. The VaR values under 95%, 97.5% and 99% confidence intervals are, respectively, explored. The results show that when the confidence interval is 95%, the investment failure rate is 4.53%; when the confidence interval is 97.5%, the investment failure rate is 2.17%; when the confidence interval is 99%, the investment failure rate is 1.08%. All these indicate that the model constructed in this study can be used in the measurement of portfolio VaR and can provide good theoretical guidance for investors.

The Shanghai Pudong Development Bank (sh600000), Guizhou Moutai (sh600519) and China Ping An Insurance (Sh601318) are taken as the research objects to explore the effect of deep learning Markov model in stock price forecasting. The number of neurons in the hidden layer (N), the coverage time of input data (t), the width of input data (d) and the training time (t) in RNN on the average prediction error of the model are compared, as shown in Figure 10. As shown in Figure 10A, with the increase in training time, the increase of input data coverage time will reduce the average prediction error of model prediction. When t = 5 and T = 16, the average prediction error is the smallest (2.89%). As shown in Figure 10B, with the increase in the number of neurons in the hidden layer of the model, the increase of input data width will also reduce the average prediction error of model prediction. When d = 5 and N = 60, the average prediction error is the smallest (2.89%). Therefore, t = 5, d = 5, T = 16 and N = 60 are selected for subsequent experiments.

Fig. 10

Subsequently, it is applied to the stock price forecast of Shanghai Pudong Development Bank (sh600000), Guizhou Moutai (sh600519) [31], and China Ping An Insurance (Sh601318), as shown in Figure 11. The variation trend of each stock price predicted by deep learning Markov model is basically consistent with the actual situation. Moreover, the errors of deep learning Markov model in stock price prediction of Shanghai Pudong Development Bank (sh600000), Guizhou Moutai (sh600519)V and China Ping An Insurance (Sh601318) are 2.56%, 2.98% [32] and 3.56%, respectively, which indicates that the deep learning Markov model proposed in this study can be used in the prediction of stock prices and has a high prediction accuracy, which is consistent with the research results of Tingwei et al. (2018) [33].

Fig. 11