SET NOECHO
PROC GAUSS
SET NOOUTPUT
* ************************************************************************
* Estimation of a function by the Gauss-Newton algorithm
*
* REFERENCE: Section 12.2.3b of Judge, Hill et al. [1988, p.514-5]
*
* SHAZAM procedure written by Diana Whistler - January 1994
*
* Input requirements:
*  (1)  NCOEF:           Number of coefficients to estimate
*  (2)  Starting values named B1, B2, .... 
*  (3)  NITER:           Maximum number of iterations
*  (4)  DELTA:           Convergence criteria - the algorithm stops when 
*                        the gradient sum of squares is less than the 
*                        value specified.
*  (5)  Y:               The left-hand side of the equation
*  (6)  EQRHS:           The right-hand side of the equation with 
*                        unknown parameters B1, B2, ...                        
*
* Run the procedure with the commands:
*   FILE PROC procedure_filename
*   EXEC GAUSS
* ************************************************************************
SET NODOECHO
GEN1 T_=$N
DIM B_ [NCOEF] BOLD_ [NCOEF] Z_ T_ [NCOEF] SSEV_ 10 STEP_ 10
DO #=1,[NCOEF]
GEN1 B_:#=B#
ENDO
DO #=1,[NITER]
* Use the DERIV command to get the Z matrix of first derivatives 
DO %=1,[NCOEF]
DERIV B% Z_:%=[EQRHS]
ENDO
* Residuals
GENR E_=[Y] - [EQRHS]
* Compute the error sum of squares
MATRIX SSE_=E_'E_
* Compute the value of the log-likelihood function
GEN1 LLF_=-(T_/2)*LOG(2*$PI)-(T_/2)-(T_/2)*LOG(SSE_/T_)
* Calculate the gradient
MATRIX G_=-2*Z_'E_
MATRIX H_=INV(Z_'Z_)
* Print results
GEN1 ITER_=#
PRINT ITER_ SSE_ LLF_ B_ G_
* Convergence test - stop if the gradient is "close" to zero.
MATRIX C_=G_'G_
ENDIF (C_.LT.[DELTA])
* Calculate the direction vector
MATRIX D_=-0.5*H_*G_
* A grid-search 
MATRIX BOLD_=B_
DO %=1,10
GEN1 STEP_:%=%/10
GEN1 STEPSZ_=STEP_:%
* Updating step
MATRIX B_=BOLD_+STEPSZ_*D_
DO !=1,[NCOEF]
GEN1 B!=B_:!
ENDO
GENR E_=[Y] - [EQRHS]
MATRIX SSE_=E_'E_
GEN1 SSEV_:%=SSE_
ENDO
* Choose an optimal step-size
SORT SSEV_ STEP_ / BEG=1 END=10
GEN1 STEPSZ_=STEP_:1
PRINT STEPSZ_
* Update the parameter vector
MATRIX B_=BOLD_+STEPSZ_*D_
DO !=1,[NCOEF]
GEN1 B!=B_:!
ENDO
ENDO
* Compute the standard errors and print the final results
GEN1 SIG2_=SSE_/T_
MATRIX SE_=SQRT(SIG2_*DIAG(H_))
MATRIX TRATIO_=B_/SE_
PRINT SSE_ LLF_ SIG2_ 
SAMPLE 1 [NCOEF]
PRINT B_ SE_ TRATIO_ 
DELETE / ALL_
SET DOECHO
SET ECHO
SET OUTPUT
PROCEND