NOTE: PROCEDURE PRINTTO used (Total process time): real time 0.00 seconds cpu time 0.00 seconds 16618 16619 proc iml; NOTE: IML Ready 16620 16621 *reset print; 16622 16623 use new.hogfuturespre1996; 16624 read all into data; 16625 16626 Num_Contracts = 6; 16627 matur = j(6,1,0); 16628 do i = 1 to 6; 16629 matur[i,] = ((i-1)*4+1)/12; 16630 end; 16631 dt = 7/360; 16632 16633 * SELECT INITIAL VALUES ; 16634 /* k = 1.489003; 16635 sigmax = 0.276;*139252; 16636 lambdax = 0.155066; 16637 mu = -0.0213; 16638 sigmae = 0.137067; 16639 rnmu = 0.0115; 16640 rho = 0.310058; 16641 psi7 = k // sigmax // lambdax // mu // sigmae // rnmu // rho;*/ 16642 s_guess = 0.0000001; 16643 16644 init_dim = 7; 16645 16646 y = data[,6:(Num_Contracts+4)]; 16646! * y is a {nobs x N} Matrix, N = number of future 16646! contracts, nobs = number of observations; 16647 nobs = nrow(y); 16648 N = ncol(y); 16649 16650 boundary = 1e300; 16650! /* bounderies for constraints */ 16651 nboundary = -1*1e300; 16652 16653 start liklhd(x0,C0) global(y,c,d,G,F,V,N,nobs,dt,psi,matur,s_guess); 16654 16655 x0 = {4.24 0}; 16655! /* Initial state vector m(t)=E[xt;et] 1x2 */ 16656 C0 = {0.1 0.1, 0.1 0.1}; 16656! /* Initial covariance matrix for the state variables 16656! W(t)=cov[xt,et];*/ 16657 m = ncol(x0); 16657! * m = Number of state variables (number of rows in a0); 16658 16659 /* THE TRANSITION EQUATION */ 16660 16661 * Extracting initial parameter values from initial psi ; 16662 16663 k = psi[1,]; 16664 sigmax = psi[2,]; 16665 lambdax = psi[3,]; 16666 mu = psi[4,]; 16667 sigmae = psi[5,]; 16668 rnmu = psi[6,]; 16669 rho = psi[7,]; 16670 16671 /* NOTATION: x(t)=c+G*x(t-1)+n(t) n~N(0,R) */ 16672 16673 c = j(m,1,0); 16673! * c is a {m x 1} Vector; 16674 c[m,] = mu*dt; 16675 G = i(m); 16675! * G is a {m x m} Matrix; 16676 G[1,1] = exp(-k*dt); 16677 16678 * Defining R = var[n(t)] and W = var[w(t)]; 16679 W11=(1-exp(-2*k*dt))*((sigmax)**2)/(2*k); 16679! /* equation 3(b) in S&S */ 16680 W12=(1-exp(-k*dt))*((rho*sigmax*sigmae)/k); 16681 W22=((sigmae)**2)*dt; 16682 W=i(m); 16683 W[1,1]=W11; 16683! W[1,2]=W12; 16683! W[2,1]=W12; 16683! W[2,2]=W22; 16684 R=i(m); 16685 16686 /* THE MEASUREMENT EQUATION */ 16687 16688 /* NOTATION: y(t)=d(t)+F(t)'x(t)+v(t) v~N(0,V) */ 16689 16690 d = j(N,1,0); 16690! * d is a {N x 1} Vector; 16691 F = j(N,m,0); 16691! * F is a {N x m} Matrix; 16692 16693 do i=1 to N; 16694 d1=rnmu*matur[i]-(1-exp(-k*matur[i]))*(lambdax/k); 16694! /* equation (9) in S&S */ 16695 d2=(1-exp(-2*k*matur[i]))*((sigmax)**2)/(2*k); 16696 d3=((sigmae)**2)*matur[i]; 16697 d4=2*(1-exp(-k*matur[i]))*((rho*sigmax*sigmae)/k); 16698 d[i,1]=d1+0.5*(d2+d3+d4); 16699 F[i,1]=exp(-k*matur[i]); 16700 F[i,2]=1; 16701 end; 16702 16703 * Measurment errors Var-Cov Matrix: Cov[v(t)]=V; 16704 16705 V_N = j(1,N,s_guess); 16706 V = diag(V_N); 16707 var = block(R,V); 16707! /* for Kalman Filter routine */ 16708 16709 nz = nrow(G); 16709! nn = nrow(y); 16709! nk = ncol(y); 16710 call kalcvf(pred,vpred,filt,vfilt,y,0,c,G,d,F,var,x0,C0); 16711 16712 et = y - pred*F`; 16713 16714 sum1 = 0; 16714! sum2 = 0; 16715 do i = 1 to nn; 16716 vpred_i = vpred[(i-1)*nz+1:i*nz,]; 16717 et_i = et[i,]; 16718 et_i1=et_i`; 16719 ft = F*vpred_i*F` + var[nz+1:nz+nk,nz+1:nz+nk]; 16720 sum1 = sum1 + log(det(ft)); 16721 sum2 = sum2 + et_i*inv(ft)*et_i1; 16722 sum = sum1 + sum2; 16723 end; 16724 return(sum); 16725 finish; NOTE: Module LIKLHD defined. 16726 16727 start loglik(y); 16728 nn = nrow(y); 16728! nn = nrow(y); 16728! nk = ncol(y); 16729 pi = constant('pi'); 16730 const = nk*nn*log(2*pi); 16731 sum = liklhd(x0,C0); 16732 log_l=(-.5*const-.5*(sum)/(nk*nn)); 16733 return(log_L); 16734 finish; NOTE: Module LOGLIK defined. 16735 16736 rank = 0; 16737 16738 do ii1 = 3.48 to 3.48 by 1; 16739 do ii2 = 0.276 to 0.336 by 0.06; 16740 do ii3 = 0.118 to 0.158 by 0.04; 16741 do ii4 = -0.0225 to 0.0125 by 0.035; 16742 do ii5 = 0.001 to 0.151 by 0.05; 16743 do ii6 = 0.00125 to 0.01125 by 0.01; 16744 do ii7 = 0.29 to 0.31 by 0.02; 16745 16746 rank = rank + 1; 16747 k = ii1; 16748 sigmax = ii2; 16749 lambdax = ii3; 16750 mu = ii4; 16751 sigmae = ii5; 16752 rnmu = ii6; 16753 rho = ii7; 16754 16755 psi7 = k // sigmax // lambdax // mu // sigmae // rnmu // rho; 16756 psi = j((init_dim+N),1,0); 16757 psi = psi7 // s_guess // s_guess // s_guess // s_guess // s_guess; 16758 16759 bounds = j(2,(init_dim+N),.); 16760 nb = ncol(bounds) - 2; 16761 lb = j((init_dim+N),1,0); 16762 lb7 = repeat(nboundary,init_dim,1); 16763 lb7[1] = 0; 16763! lb7[2] = 0; 16763! lb7[5] = 0; 16763! lb7[7] = -1; 16764 lb5 = repeat(0.000000001,5,1); 16765 lb = lb7 // lb5; 16766 16767 ub = j((init_dim+N),1,0); 16768 ub6 = repeat(boundary,6,1); 16769 ub7 = 1; 16770 ub5 = repeat(boundary,5,1); 16771 ub = ub6 // ub7 // ub5; 16772 do jj = 1 to nb; 16773 bounds[1,jj] = lb[jj,]; 16774 bounds[2,jj] = ub[jj,]; 16775 end; 16776 16777 16778 opt = {0 5 . 1}; 16779 tc = {2000 5000}; 16780 call nlpqn(rc,psi_opt, "loglik", psi, opt, bounds, tc,,,,,); 16781 psi_opt1 = t(psi_opt); 16782 sumd = 0; 16783 do i = 1 to 7; 16784 sumd = sumd + abs(psi[i,] - psi_opt[i]); 16785 end; 16786 /*if sumd > 0 then 16786! */do; 16787 print rank sumd psi_opt1 psi pred; 16788 end; 16789 16790 end; 16791 end; 16792 end; 16793 end; 16794 end; 16795 end; 16796 end; NOTE: ABSGCONV convergence criterion satisfied. ERROR: Matrix pred has not been set to a value. statement : PRINT at line 16787 column 3 16797 16798 quit; NOTE: Exiting IML. NOTE: The SAS System stopped processing this step because of errors. NOTE: PROCEDURE IML used (Total process time): real time 0.67 seconds cpu time 0.34 seconds 16799 16800 proc printto; 16801 run;