'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