我想做一个BEKK模型，原来已经弄好了程序，后来又想把正态分布假设改成 t 分布，但是在编程上t分布就不会了，希望有哪位好心人能帮帮我
原程序如下：
smpl @all
series y1 =r1
series y2 =r2

sample s0 1 692
sample s1 2 692

smpl s0
equation eq1.arch(m=100,c=1e-5) y1  c
equation eq2.arch(m=100,c=1e-5) y2 c
coef(2) mu
mu(1) = eq1.c(1)
mu(2)= eq2.c(1)

coef(3) omega
omega(1)=(eq1.c(2))^.5
omega(2)=0
omega(3)=eq2.c(2)^.5

coef(4) alpha
alpha(1) = (eq1.c(3))^.5
alpha(2) = (eq2.c(3))^.5
alpha(3) = 0
alpha(4) = 0

coef(5) beta
beta(1) = (eq1.c(4))^.5
beta(2) = (eq2.c(4))^.5
beta(3) = 0
beta(4) = 0

!mlog2pi = 2*log(2*@acos(-1))

series cov_y1y2 = @cov(y1-mu(1), y2-mu(2))
series var_y1 = @var(y1)
series var_y2 = @var(y2)

series sqres1 = (y1-mu(1))^2
series sqres2 = (y2-mu(2))^2
series res1res2 = (y1-mu(1))*(y2-mu(2))  

logl bvgarch
bvgarch.append @logl logl
bvgarch.append sqres1 = (y1-mu(1))^2
bvgarch.append sqres2 = (y2-mu(2))^2
bvgarch.append res1res2 = (y1-mu(1))*(y2-mu(2))

bvgarch.append var_y1  =  omega(1)^2 + beta(1)^2*var_y1(-1) + 2*beta(1)*beta(4)*cov_y1y2(-1) + beta(4)^2*var_y2(-1)+ alpha(1)^2*sqres1(-1) + 2*alpha(1)*alpha(4)*res1res2(-1) + alpha(4)^2*sqres2(-1)

bvgarch.append var_y2  = omega(3)^2+omega(2)^2 + beta(3)^2*var_y1(-1) + 2*beta(2)*beta(3)*cov_y1y2(-1) + beta(2)^2*var_y2(-1)+ alpha(3)^2*sqres1(-1) + 2*alpha(2)*alpha(3)*res1res2(-1) + alpha(2)^2*sqres2(-1)

bvgarch.append cov_y1y2 = omega(1)*omega(2) + beta(1)*beta(3)*var_y1(-1) + beta(2)*beta(4)*var_y2(-1) + beta(3)*beta(4)*cov_y1y2(-1)  + beta(1)*beta(2)*cov_y1y2(-1)  + alpha(1)*alpha(3)*sqres1(-1) + alpha(3)*alpha(4)*res1res2(-1) +alpha(1)*alpha(2)*res1res2(-1) + alpha(2)*alpha(4)*sqres2(-1)

bvgarch.append deth = var_y1*var_y2 - cov_y1y2^2

' inverse elements of the variance-covariance matrix
bvgarch.append invh1 = var_y2/deth
bvgarch.append invh3 = var_y1/deth
bvgarch.append invh2 = -cov_y1y2/deth

' log-likelihood series
bvgarch.append logl =-0.5*(!mlog2pi + (invh1*sqres1+2*invh2*res1res2+invh3*sqres2) + log(deth))

' remove some of the intermediary series
' bvgarch.append @temp invh1 invh2 invh3 sqres1 sqres2 res1res2 deth

smpl s1
bvgarch.ml(showopts, m=100, c=1e-5)

show bvgarch.output
delete y1 y2 alpha beta omega deth mu invh1 invh2 invh3 logl cov_y1y2 res1res2 var_y1 var_y2 sqres1 sqres2 s0 s1 eq1 eq2