home *** CD-ROM | disk | FTP | other *** search
/ Source Code 1992 March / Source_Code_CD-ROM_Walnut_Creek_March_1992.iso / usenet / altsrcs / 2 / 2400 / nlopr.cc next >
Text File  |  1990-12-28  |  9KB  |  525 lines
  1. /* ---------------------------------------------------------------------------
  2.  
  3. nlmdl: nlopr.cc
  4.  
  5. nlmdl is a C++ implementation of the statistical methods in A. Ronald 
  6. Gallant, "Nonlinear Statistical Models," New York: John Wiley and Sons, 
  7. 1987, ISBN 0-471-80260-3, using a matrix class realmat that is distributed 
  8. with it.  The header files nlmdl.h and realmat.h describe the use of the 
  9. program and matrix class, respectively.  
  10.  
  11. Copyright (C) 1990.
  12.  
  13. A. Ronald Gallant
  14. P.O. Box 5513 
  15. Raleigh NC 27650-5513 
  16. USA   
  17.  
  18. Permission to use, copy, modify, and distribute this software and its 
  19. documentation for any purpose and without fee is hereby granted, provided 
  20. that the above copyright notice appear in all copies and that both that 
  21. copyright notice and this permission notice appear in supporting 
  22. documentation.  
  23.  
  24. This software is provided "as is" without any expressed or implied warranty.
  25.  
  26. ------------------------------------------------------------------------------
  27.  
  28. This is the collection of operators used by nlmdl.  They act on s, of class 
  29. status, using m, of class model; s must have been filled in and the function 
  30. initialize(), of class model, must have been called before using any of these 
  31. operators.  
  32.  
  33. --------------------------------------------------------------------------- */
  34.  
  35. int line_search(char* msg)
  36. {
  37.   REAL    norm_D, norm_theta;
  38.   REAL    step_length;
  39.   REAL    save_obj = s.obj;
  40.   realmat save_theta = s.theta;
  41.   realmat ss(1,1);
  42.   
  43.   ss = T(s.D) * s.D;
  44.   norm_D = sqrt(ss[1]);
  45.  
  46.   ss = T(s.theta) * s.theta;
  47.   norm_theta = sqrt(ss[1]);
  48.  
  49.   if ( norm_D < s.tol*(norm_theta + s.eps) ) {
  50.     strcpy(msg,"Tolerence check passed.\n");
  51.     return 1;
  52.   }
  53.  
  54.   for (step_length=1; step_length >= 0.5; step_length -= 0.1) {
  55.  
  56.     s.theta = save_theta + step_length * s.D;
  57.  
  58.     (*opr_obj)();
  59.  
  60.     if (s.obj < save_obj) {
  61.       sprintf(msg,"Step length = %g \n",step_length);
  62.       return 0 ;
  63.     }
  64.   }
  65.  
  66.   for (step_length=0.5; step_length > 1.0e-13; step_length *= 0.5) {
  67.  
  68.     s.theta = save_theta + step_length * s.D;
  69.  
  70.     (*opr_obj)();
  71.  
  72.     if (s.obj < save_obj) {
  73.       sprintf(msg,"Step length = %g \n",step_length);
  74.       return 0;
  75.     }
  76.   }
  77.  
  78.   strcpy(msg,"A line search did not improve this estimate.\n");
  79.   s.obj = save_obj;
  80.   s.theta = save_theta;
  81.   return 1;
  82. }
  83.  
  84.  
  85. void SUR_obj()
  86. {
  87.   INTEGER t; 
  88.   realmat sse(1,1,(REAL)0 );
  89.   realmat varinv = invpsd(s.var,s.eps);
  90.  
  91.   for (t=1; t<=s.n; t++) 
  92.     sse = sse + T(m.e(t)) * varinv * m.e(t);
  93.  
  94.   s.obj = sse[1];
  95. }
  96.  
  97.  
  98. void SUR_mgn()
  99. {
  100.   INTEGER t; 
  101.   realmat sse;
  102.   realmat varinv = invpsd(s.var,s.eps);
  103.  
  104.   sse.resize(1,1,(REAL)0 );
  105.   s.V.resize(s.p, s.p, (REAL)0 );
  106.   s.D.resize(s.p, 1, (REAL)0 );
  107.  
  108.   for(t=1; t<=s.n; t++) {
  109.     sse = sse + T(m.e(t))    * varinv * m.e(t);
  110.     s.V = s.V + T(m.dele(t)) * varinv * m.dele(t);
  111.     s.D = s.D + T(m.dele(t)) * varinv * m.e(t);
  112.   }
  113.  
  114.   s.obj = sse[1];
  115.   s.V = invpsd(s.V,s.eps);
  116.   s.D = - s.V * s.D;
  117. }
  118.  
  119.  
  120. void SUR_var(int var_loop)
  121. {
  122.   INTEGER i,j,t;
  123.   REAL fn;
  124.   
  125.   if (rows(s.var)==0 || cols(s.var)==0) {
  126.     s.var.resize(s.M,s.M,(REAL)0 );
  127.     for (i=1; i<=s.M; i++) s.var.elem(i,i)=1.0;
  128.   }
  129.  
  130.   if (var_loop == 0) return;
  131.  
  132.   for (i=1; i<=s.M; i++) 
  133.   for (j=1; j<=s.M; j++) 
  134.     s.var.elem(i,j) = (REAL)0 ;
  135.  
  136.   for (t=1; t<=s.n; t++) 
  137.     s.var = s.var + m.e(t) * T(m.e(t));
  138.  
  139.   if (s.M == 1 && s.n > s.p) {
  140.     fn = s.n - s.p;
  141.     strcpy(s.df,"corrected");
  142.   }
  143.   else {
  144.     fn = s.n;
  145.     strcpy(s.df,"uncorrected");
  146.   }
  147.  
  148.   for (i=1; i<=s.M; i++) 
  149.   for (j=1; j<=s.M; j++) 
  150.     s.var.elem(i,j) = s.var.elem(i,j)/fn;
  151. }
  152.  
  153.  
  154. void SUR_V() 
  155. { 
  156.   INTEGER t; 
  157.   INTEGER tau; 
  158.   INTEGER i; 
  159.   REAL    x,weight;
  160.   realmat varinv = invpsd(s.var,s.eps);
  161.   
  162.   if (strcmp(s.vartype,"heteroskedastic") == 0 || s.MA > 0) {
  163.  
  164.     realmat I(s.p,s.p,(REAL)0 );
  165.     realmat I_tau;
  166.  
  167.     for (tau=0; tau<=s.MA; tau++) {
  168.  
  169.       I_tau.resize(s.p,s.p,(REAL)0 );
  170.  
  171.       for (t=tau+1; t<=s.n; t++)
  172.         I_tau = I_tau 
  173.                 + T(m.dele(t)) * varinv * m.e(t)  
  174.                   * T(m.e(t-tau)) * varinv * m.dele(t-tau);
  175.  
  176.       if (strcmp(s.weights,"Parzen") == 0 && tau > 0) {
  177.         x = tau/(REAL)s.MA;
  178.  
  179.         if ( x < 0.5 ) 
  180.           weight = 1.0 - 6.0*pow(x,2) + 6.0*pow(x,3);
  181.         else
  182.           weight = 2.0*pow((1.0 - x),3);
  183.       }
  184.       else {
  185.         weight = 1.0;
  186.       }
  187.  
  188.       I = I + weight*I_tau;
  189.  
  190.       if (tau > 0)
  191.         I = I + weight*T(I_tau);
  192.     }
  193.  
  194.     s.V = s.V * I * s.V;
  195.   }
  196.  
  197.   s.rank = 0;
  198.   for (i=1; i<=s.p; i++)
  199.     if (s.V.elem(i,i) > (REAL)0  ) s.rank++;
  200.  
  201.   return;
  202. }
  203.  
  204.  
  205. realmat qZ;
  206.  
  207. void opr_qZ(INTEGER t)
  208.  
  209. {
  210.   INTEGER alpha, i, ii;
  211.  
  212.   realmat q_tmp(s.M,1);
  213.   realmat Z_tmp(s.K,1);
  214.  
  215.   qZ.resize(s.M*s.K,1);
  216.  
  217.   q_tmp = m.e(t);
  218.   Z_tmp = m.Z(t);
  219.  
  220.   for (alpha=1; alpha<=s.M; alpha++) {
  221.   for (i=1; i<=s.K; i++)             {
  222.  
  223.     ii = s.K*(alpha-1) + i;
  224.  
  225.     qZ[ii] = q_tmp[alpha] * Z_tmp[i];
  226.   }
  227.   }
  228. }
  229.  
  230.  
  231. realmat QZ;
  232.  
  233. void opr_QZ(INTEGER t)
  234. {
  235.   INTEGER alpha, i, ii, k;
  236.  
  237.   realmat Q_tmp(s.M,s.p);
  238.   realmat Z_tmp(s.K,1);
  239.  
  240.   QZ.resize(s.M*s.K,s.p);
  241.  
  242.   Q_tmp = m.dele(t);
  243.   Z_tmp = m.Z(t);
  244.  
  245.   for (k=1; k<=s.p; k++)             {
  246.   for (alpha=1; alpha<=s.M; alpha++) {
  247.   for (i=1; i<=s.K; i++)             {
  248.  
  249.     ii = s.K*(alpha-1) + i;
  250.  
  251.     QZ.elem(ii,k) = Q_tmp.elem(alpha,k) * Z_tmp[i];
  252.   }
  253.   }
  254.   }
  255. }
  256.  
  257.  
  258. realmat qZ_mts;
  259.  
  260. void opr_qZ_mts()
  261. { 
  262.  
  263.   INTEGER t;
  264.   
  265.   qZ_mts.resize(s.M*s.K, 1, (REAL)0 );
  266.  
  267.   for (t=1; t<=s.n; t++) {
  268.  
  269.     opr_qZ(t);
  270.  
  271.     qZ_mts = qZ_mts + qZ;
  272.  
  273.   }
  274. }
  275.  
  276.  
  277. realmat QZ_mts;
  278.  
  279. void opr_QZ_mts()
  280. { 
  281.   INTEGER t;
  282.   
  283.   QZ_mts.resize(s.M*s.K, s.p, (REAL)0 );
  284.  
  285.   for (t=1; t<=s.n; t++) {
  286.  
  287.     opr_QZ(t);
  288.  
  289.     QZ_mts = QZ_mts + QZ;
  290.  
  291.   }
  292. }
  293.  
  294.  
  295. realmat ZZ;
  296.  
  297. void opr_ZZ()
  298. {
  299.   INTEGER t;
  300.  
  301.   ZZ.resize(s.K, s.K, (REAL)0 );
  302.  
  303.   for (t=1; t<=s.n; t++)
  304.     ZZ = ZZ + m.Z(t) * T(m.Z(t));
  305. }
  306.  
  307.  
  308. void TSLS_obj() 
  309. {
  310.  
  311.   INTEGER alpha, i, ii;
  312.   INTEGER beta,  j, jj;
  313.  
  314.   realmat varinv = invpsd(s.var,s.eps);
  315.  
  316.   opr_qZ_mts();
  317.   opr_ZZ();
  318.  
  319.   realmat ZZinv = invpsd(ZZ,s.eps);
  320.   s.obj = (REAL)0 ;
  321.  
  322.   for (alpha=1; alpha<=s.M; alpha++) {
  323.   for (i=1; i<=s.K; i++)             {
  324.   for (beta=1; beta<=s.M; beta++)    {
  325.   for (j=1; j<=s.K; j++)             {
  326.  
  327.     ii = s.K*(alpha-1) + i;
  328.     jj = s.K*(beta-1)  + j;
  329.  
  330.     s.obj += qZ_mts[ii] * qZ_mts[jj] 
  331.              * varinv.elem(alpha,beta) * ZZinv.elem(i,j);
  332.   }
  333.   }
  334.   }
  335.   }
  336. }
  337.  
  338.  
  339. void TSLS_mgn() 
  340. { 
  341.   INTEGER alpha, i, ii;
  342.   INTEGER beta,  j, jj;
  343.   INTEGER k,l;
  344.  
  345.   realmat varinv = invpsd(s.var,s.eps);
  346.  
  347.   opr_qZ_mts();
  348.   opr_QZ_mts();
  349.   opr_ZZ();
  350.  
  351.   realmat ZZinv = invpsd(ZZ,s.eps);
  352.  
  353.   s.obj = (REAL)0 ;
  354.   s.V.resize(s.p, s.p, (REAL)0 );
  355.   s.D.resize(s.p, 1, (REAL)0 );
  356.  
  357.   for (alpha=1; alpha<=s.M; alpha++) {
  358.   for (i=1; i<=s.K; i++)             {
  359.   for (beta=1; beta<=s.M; beta++)    {
  360.   for (j=1; j<=s.K; j++)             {
  361.  
  362.     ii = s.K*(alpha-1) + i;
  363.     jj = s.K*(beta-1)  + j;
  364.  
  365.     s.obj += qZ_mts[ii] * qZ_mts[jj] 
  366.              * varinv.elem(alpha,beta) * ZZinv.elem(i,j);
  367.  
  368.   for (k=1; k<=s.p; k++)             {
  369.  
  370.     s.D[k] += QZ_mts.elem(ii,k) * qZ_mts[jj] 
  371.               * varinv.elem(alpha,beta) * ZZinv.elem(i,j);
  372.  
  373.   for (l=1; l<=s.p; l++)             {
  374.  
  375.     s.V.elem(k,l) += QZ_mts.elem(ii,k) * QZ_mts.elem(jj,l) 
  376.                      * varinv.elem(alpha,beta) * ZZinv.elem(i,j);
  377.  
  378.   }
  379.   }
  380.   }
  381.   }
  382.   }
  383.   }
  384.  
  385.   s.V = invpsd(s.V,s.eps);
  386.  
  387.   s.D = - s.V * s.D;
  388.  
  389. }
  390.  
  391.  
  392. void TSLS_V() 
  393. { 
  394.   INTEGER i;
  395.  
  396.   s.rank = 0;
  397.   for (i=1; i<=s.p; i++)
  398.     if (s.V.elem(i,i) > (REAL)0  ) s.rank++;
  399.  
  400.   return;
  401. };
  402.  
  403.  
  404. void GMM_obj() 
  405. { 
  406.  
  407.   realmat sse(1,1);
  408.   realmat varinv = invpsd(s.var,s.eps);
  409.  
  410.   opr_qZ_mts();
  411.  
  412.   sse = T(qZ_mts) * varinv * qZ_mts ;
  413.  
  414.   s.obj = sse[1];
  415. }
  416.  
  417.  
  418. void GMM_mgn() { 
  419.  
  420.   realmat sse(1,1);
  421.   realmat varinv = invpsd(s.var,s.eps);
  422.  
  423.   opr_qZ_mts();
  424.   opr_QZ_mts();
  425.  
  426.   sse = T(qZ_mts) * varinv * qZ_mts ;
  427.  
  428.   s.obj = sse[1];
  429.  
  430.   s.D = T(QZ_mts) * varinv * qZ_mts ;
  431.  
  432.   s.V = T(QZ_mts) * varinv * QZ_mts ;
  433.  
  434.   s.V = invpsd(s.V,s.eps);
  435.  
  436.   s.D = - s.V * s.D;
  437.  
  438. }
  439.  
  440.  
  441. void GMM_var(int var_loop) 
  442. {
  443.   INTEGER l = s.M * s.K;
  444.   INTEGER t; 
  445.   INTEGER tau; 
  446.   INTEGER alpha, i, j, ii, jj; 
  447.   REAL    x,weight;
  448.  
  449.   strcpy(s.df,"uncorrected");
  450.  
  451.   if (rows(s.var)==0 || cols(s.var)==0) {
  452.  
  453.     s.var.resize(l,l,(REAL)0 );
  454.  
  455.     opr_ZZ();
  456.  
  457.     for (alpha=1; alpha<=s.M; alpha++) {
  458.     for (i=1; i<=s.K; i++)             {
  459.     for (j=1; j<=s.K; j++)             {
  460.  
  461.       ii = s.K*(alpha-1) + i;
  462.       jj = s.K*(alpha-1) + j;
  463.  
  464.       s.var.elem(ii,jj) = ZZ.elem(i,j);
  465.  
  466.     }
  467.     }
  468.     }
  469.   }
  470.   if (var_loop == 0) return;
  471.  
  472.   realmat I(l,l);
  473.   realmat I_tau(l,l);
  474.   realmat qZ_lag(l,1);
  475.  
  476.   for (tau=0; tau<=s.MA; tau++) {
  477.  
  478.     for (i=1; i<=l; i++)
  479.     for (j=1; j<=l; j++)
  480.       I_tau.elem(i,j)=(REAL)0 ;
  481.  
  482.     for (t=tau+1; t<=s.n; t++) {
  483.  
  484.       opr_qZ(t-tau);
  485.       qZ_lag = qZ;
  486.       opr_qZ(t);
  487.       
  488.       I_tau = I_tau + qZ * T(qZ_lag);
  489.     }
  490.  
  491.     if (strcmp(s.weights,"Parzen") == 0 && tau > 0) {
  492.       x = tau/(REAL)s.MA;
  493.       if ( x < 0.5 ) 
  494.         weight = 1.0 - 6.0*pow(x,2) + 6.0*pow(x,3);
  495.       else
  496.         weight = 2.0*pow((1.0 - x),3);
  497.         
  498.     }
  499.     else {
  500.       weight = 1.0;
  501.     }
  502.  
  503.     I = I + weight*I_tau;
  504.     if (tau > 0) {
  505.       I = I + weight*T(I_tau);
  506.     }
  507.   }
  508.  
  509.   s.var = I;
  510.  
  511.   return;
  512. }
  513.  
  514.  
  515. void GMM_V() 
  516. { 
  517.   INTEGER i;
  518.  
  519.   s.rank = 0;
  520.   for (i=1; i<=s.p; i++)
  521.     if (s.V.elem(i,i) > (REAL)0  ) s.rank++;
  522.  
  523.   return;
  524. }
  525.