'Two independent variable regression model with t-distributed errors

'Create a workfile
wfcreate ytmodel u 1000

'Define the coefficients
!rho1 = 0.6
!rho2 = 0.3
!dof = 3

'Generate the series
rndseed 12345
series x1 = @rnorm
series x2 = @rnorm
series y = !rho1*x1 + !rho2*x2 + @rtdist(!dof)

'Declare coefficent vectors and assign initial values
coef(2) rho = .5
coef(1) dof = 2
coef(1) var = 1
!pi = @acos(-1)

'Set up logl object
logl yt
yt.append @logl maxlik
yt.append res = y - rho(1)*x1 - rho(2)*x2
yt.append sqres = ((res^2/var(1))/dof(1))+ 1
yt.append maxlik = @gammalog((dof(1) + 1)/2) - @gammalog(dof(1)/2) - log(!pi)/2 - log(dof(1))/2 - log(var(1))/2 - (dof(1)+1)*log(sqres)/2

'Run MLE and display results
smpl @all
yt.ml(showopts,m=1000,c=1e-5)
show yt.output

'Joint test for parameter values
yt.wald rho(1)=!rho1,rho(2)=!rho2,dof(1)=!dof

'Put the model into state space form and estimate in Gaussian form
sspace ss_normal
ss_normal.append @signal y = sv1*x1 + sv2*x2 + [var=exp(c(1))]
ss_normal.append @state sv1 = sv1(-1)
ss_normal.append @state sv2 = sv2(-1)
ss_normal.append param c(1) .0
ss_normal.ml
show ss_normal

'Estimate the model via add-in
ss_normal.sspacetdist(iters=50,delta=1e-02,adjsave)
show ss_normal_new

'Distribution fit for disturbance term
ss_normal_new.makesignals(t=disturb) disturbs
freeze(mode=overwrite,distfit) disturbs.distplot hist(anchor=0, scale=dens) theory(lgnd=detail) theory(dist=tdist, lgnd=detail)
show distfit