MATRIX - Section 52 The B34S programming language capability is called with the command b34sexec matrix; /$ commands here b34srun; Example: b34sexec matrix; x=matrix(2,2:22. 33. 44. .02); ix=inv(x); test=ix*x; call print(x,ix,test); b34srun; The MATRIX command provides a 4th generation language to process selected calculations from B34S procedures and the ability to further process data. Analysis is supported for real*8, real*16, complex*16 and complex*32 data. Other data types such as character*1, character*8, real*4 and integer*4 are supported. Calculation of the inverse of a real*4 matrix, while supported, is not recommend due to accuracy loss. High resolution graphics are available on all currently supported platforms (Windows, RS/6000, Sun, Linux) and batch and interactive operation is supported. Many string operations are available for character*8 and character*1 data. The MATRIX facility supports user PROGRAMS, SUBROUTINES and FUNCTIONS. SUBROUTINES and FUNCTIONS have their own address space. Variables are built using object oriented programming with analytic statements such as: y = x*2.; where x is a variable that could be a matrix, 2D array, 1D array, vector or a scalar. The use of the MATRIX command FORMULA and the SOLVE subroutine allows recursive solution of an analytic statement over a range of index values. This facility speeds up calculations that would have had to use do loops which have substantial overhead. For intensely recursive calculations, the user can call Fortran routines from within a MATRIX command program or subroutine. The MATRIX command recognizes variables of KIND complex*16, real*8, real*4, real*16, complex*32, integer*4, character*1 and character*8. Due to possible accuray loss only array math is supported for integer*4 objects. While the inv(real4object) is supported, due to accuracy loss this is not recommended. The KLASS of an object determines how it is processed. KLASS types are scalar, 1D array, 2D array, vector and 2D matrix. For example the assingment statement y=2.0; sets y to be a real*8 scalar, while y=2; sets y to be a integer*4 scalar. If the commands x=33.; y=2; test=x*y; are given there will be a mixed-mode error message since the MATRIX command processor does not know whether test should be integer*4 or real*8. To create an integer*4 from a real*8 use testint=idint(x)*y; while to create a real*8 from an integer*4 variable y use testr8 =x*dfloat(y); Where possible Fortran names have been used for built-in functions to reduce the learning curve. To create a 2D array use xarray=array(2,2: 1., 2., 3., 4.); To create a 2D matrix use xmatrix=matrix(2,2:1., 2., 3., 4.); or xmatrix=mfam(xarray); There are a number of built in functions to convert the KIND and KLASS. Key ones are: SFAM( ) Create a scalar family. MFAM( ) Convert an array to a matrix. AFAM( ) Convert a matrix to an array. VFAM( ) Convert a 1D array to a 1D vector. TO_RMATRIX( ) Convert Object to Row-Matrix TO_CMATRIX( ) Convert Object to Col-Matrix TO_RARRAY( ) Convert Object to Row-Array TO_CARRAY( ) Convert Object to Col-Matrix TO_VECTOR( ) Convert Object to Vector TO_ARRAY Convert Object to Array DFLOAT( ) Convert integer to real*8. IDINT( ) Convert real*8 to integer. IQINT( ) Convert real*16 to integer. IDNINT( ) Nearest integer from real*8. IQNINT( ) Nearest integer from real*16 REAL( ) Obtain real*8 part of a complex*16 data type. QREAL( ) Obtain real*16 part of complex*32 data type. IMAG( ) Obtain the real*8 imaginary part of a complex*16 data type. QIMAG( ) Obtain the real*16 imaginary part of a complex*32 data type. COMPLEX( , ) Build a complex*16 variable from two real*8 variables. QCOMPLEX( , ) Build a complex*32 variable from two real*16 variables. QNINT( ) Real*16 representation of nearest whole number. DNINT( ) Real*8 representation of nearest whole number. R8TOR16( ) Real*8 to real*16. R16TOR8( ) Real*16 to real*8. C16TOC32( ) Complex*16 to complex*32. C32TOC16( ) Complex*32 to complex*16. KINDAS(x,yy) Set a variable the kind of x, but value of yy. More detail on these features will be given below. The B34S MATRIX language is very close to the Speakeasy language. In MATLAB all 2D and 1D objects are assumed to be what B34S calls the MATRIX KLASS. Assume x and y are 3 by 3 MATRIX objects in B34S and MATLAB. In MATLAB for an element by element operation the command is newobj=x.*y; In B34S and Speakeasy the command for the same operation is: newobj=afam(x)*afam(y); If it is desired to place this back in a matrix: newobj=mfam(afam(x)*afam(y)); Unlike MATLAB, B34S and Speakeasy are not case sensitive. Assume xarray is an array and xmatrix is a matrix. xxarray =2.+xarray; xxmatrix=2.+xmatrix; Speakeasy and B34S work the same. Here 2. is added to all elements in xarray and only the diagonal elements in xmatrix. In MATLAB 2+x adds 2 to all elements of the matrix. The development of the MATRIX facility has been influenced closely by Speakeasy. Wherever possible, the same command language has been used to facilitate transfer of data to Speakeasy for further processing. The MATRIX facility is NOT designed to run interactively but commands can be given interactively in the MANUAL mode. In this mode scripts can be edited and submitted. Output is written to the b34s.out file and error messages are displayed in both the b34s.out and b34s.log files. Under the Display Manager, the user can scroll the output files and run the MATRIX command before and after other b34s commands. The objective of the MATRIX facility is to give the user access to a powerful object oriented programming language so that custom calculations can be made. Of major interest is providing the ability to estimate complex nonlinear least squares and maximum likelihood models. Such models, which are specified in B34S MATRIX command programs, can be solved using either other B34S subroutines or with the MATRIX nonlinear commands NLLSQ, NL2SOL, MAXF1, MAXF2, MAXF3, CMAXF1, CMAXF2 and CMAXF3. While the use of B34S subroutines would give the user TOTAL control of the estimation process, speed would be given up. In addition to specification of the model, in the B34S MATRIX language it is also possible to write the model in a Fortran or C program and call this user program from within a B34S MATRIX program. For recursive systems where it is near impossible to vectorize the calculation, this is may be the best way to proceed. The B34S nonlinear solvers are based on time tested routines. The objective of the b34s implementation is to facilitate their use in a way whereby there is full knowledge of just what is being calculated. The design of the MATRIX facility allows other libraries of routines to be accessed from C or Fortran to provide other alternatives. The MATRIX nonlinear commands give the user complete control of the form of the estimated model which is specified in a MATRIX command PROGRAM. Since these programs are called by compiled solvers, there is a substantial speed advantage over a design that writes the solver in a subroutine. The file MATRIX2.MAC contains the subroutines DUD and MARQ which illustrate the subroutine approach to nonlinear least squares and the power of the MATRIX command language. While these subroutines can be used, the NLLSQ and NL2SOL commands are substantially more powerful and orders of magnitide faster. The MATRIX command display routines OUTSTRING, OUTINTEGER and OUTDOUBLE can be used to "instrument" the solution process such that the user can see how the search is proceeding if B34S is running in the Display Manager. These commands will not work for batch jobs since the proper windows have not been opened unless the b34s2 procedure is used. The B34S Matrix Language design allows DO loops and IF structures to be in the command stream provided the commands are not given in manual mode. This is not possible with Speakeasy. The command call manual; can be placed anywhere in the job to place the processor in manual mode whereby commands could be entered or the process modified. The main purpose of manual mode is to allow the user to take control of how a subroutine is running. If call manual; is used in a MATRIX program under the Display Manager, the user is placed in the IMATRIX mode where the output to date can be seen, errors can be trapped and commands can be given to see what has been calculated to date and scripts can be submitted. If the user is debugging a subroutine it is often useful to be able to see what variables are active and possible modify the course of the subroutines execution in a real-time mode. The command call break('Are at position a'); can be used to trap execution in a loop. If any key has been hit prior to call break being found, the program will stop. Otherwise the command is ignored. This saves having the program checking for a key stroke as each command is being executed which will slow speed. If call break is executed, the user can then either allow execution to proceed or terminate the process then and there. If the MATRIX workspace is saved with the SAVE command, it can be restored after more B34S commands have run. Unlike most other B34S commands, in many cases the MATRIX command requires that commas be used. Unlike other B34S commands, errors are written both in the LOG file and at the exact command that has caused a problem. The command call echooff; can be used to reduce output when running a user SUBROUTINE, FUNCTION or PROGRAM. If there is some question on how a certain section is running, the command call echoon; can be used to echo commands as they are executed. In Matlab commands that end with ; are not echoed, while commands that do not end in ; are. This "design" requires the user to modify code during the debug phase. As a result errors can creep in. The echoon and echoff commands allow one tyo globally turn off output. Every attempt has been made to increase the speed of execution. Unlike Speakeasy, the CALL command is needed to execute PROGRAMS or SUBROUTINES. By this design, a great deal of search time is saved. To avoid I/O and to avoid possible unforseen conflicts over names, user SUBROUTINES, PROGRAMS or FUNCTIONS must be loaded to be found. Automatic loading of user commands slows execution substantially since the processor would have to look in all libraries to see if a built-in command had been replaced by a user subroutine, program or function. A major reason automatic loading was not implemented was to avoid the danger that the wrong subroutine could be loaded and a "wrong" calculation be made and not caught. A major objective of the MATRIX command is to give the user total control over the processing of data. Users can develope extensions and modifications of the statistics calculated by other B34S commands. Over the years, more and more of the regular commands will be converted to save data for later use in the matrix command. The library of applications SUBROUTINES, FUNCTIONS and PROGRAMS in c:\b34slm\matrix2.mac will be inhanced. Since these procedures are written in the B34S MATRIX language, they are self documenting. Facilities are provided whereby users can add to these libraries. Since all procedures have to be loaded, there is NO possibility that a user procedure can conflict with a built-in command unless it is explicitly loaded. The only conflict possible is with a MATRIX language command. This design prevents the MATLAB problem of suddenly key commands not working because the order of the library is changed due to a new toolkit being added. The routines in matrix2.mac are of two types. General purpose routines, such CFREQ which will calculate a cumulative frequency distribution, are documented BOTH in the TOOLKIT section which lists all routines and under the matrix command subroutine list. Commands such as this have test problems of the same name in matrix.mac. Other programs in matrix2.mac of less general interest for production jobs are only documented in the TOOLKIT section. The library staging.mac contains example files for the proposed subroutines in staging2.mac. The MATRIX command can read and process a binary file. Reading can be done in any order. This facility allows the system to process and change a load module or recover data from a file that has a strange structure. Character*1 processing allows the user to read and parse lines of a file. This capability is of interest to the expert user. For NLLSQ and maximum likelihood models that involve recursive models, in many cases an alternative program may be RATS. Such models are slow to estimate since there is a heavy parse overhead and the the vector capability in the B34S MATRIX facility is not usuable. B34S is not designed to replace RATS in this area. In B34S NLLSQ and maximum likelihood estimation can be done where the model is specified in the user's program using MATRIX commands. This mode of operation combines the advantages of fast execution with the flexibility of allowing the user to specify the model using a 4th generation language. NLLSQ can also be done using subroutines DUD and MARQ which were written in the B34S MATRIX command language. The estimation of a ML model when a recursive system is not required is fast and a number of routines from IMSL are supplied. Supported capability includes constrained ML models. While recursive systems are possible, due to the fact that DO loops or SOLVE/FORMULA commands are needed, the speed is slower. The subroutine GARCH can be used to estimate a subset of the ARCH/GARCH class of models without some of the recursive call overhead. An even simplier command GARCHEST allows univariate estimation for simple ARCH/GARCH models. Since B34S can estimate a nonlinear programming model with nonlinear constraints, a fairly wide class of models can be estimated. The BGARCH subroutine allows the user's PROGRAM to estimate a bivariate GARCH model without the recursive overhead. Note: In the above NLLSQ refers to nonlinear least squares. The B34S MATRIX command has two commands for this calculation. NLLSQ uses the Meeter (1964a, 1964b) routines that were first developed at the University of Wisconsin and form the basis of much work in time series analysis. The alternative routine NL2SOL uses the Dennis-Gay-Welsch (1981) programs that have been widely used in the literature. As of May 2004, the NLLSQ command can be run with real*8 or real*16 data. The maximization routines supported all come from various versions of the IMSL library. The B34S MATRIX command subroutines DUD and MARQ use the logic of the SAS Nonlin command but are written in the B34S command language. These are supplied for research interest. For nonlinear least squares, these would be a third choice. ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ List of Built-In Matrix Command Subroutines ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ACEFIT - Alternating Conditional Expectation Model Estimation ACF_PLOT - Simple ACF Plot ADDCOL - Add a column to a 2d array or matrix. ADDROW - Add a row to a 2d array or matrix. AGGDATA - Aggregate Data under control of an ID Vector. ALIGN - Align Series with Missing Data ARMA - ARMA estimation using ML and MOM. AUTOBJ - Automatic Estimation of Box-Jenkins Model BACKSPACE - Backspace a unit BDS - BDS Nonlinearity test. BESTREG - Best OLS REGRESSION B_G_TEST - Breusch-Godfrey (1978) Residual Test BGARCH - Calculate function for a BGARCH model. BLUS - BLUS Residual Analysis BPFILTER - Baxter-King Filter. BREAK - Set User Program Break Point. BUILDLAG - Builds NEWY and NEWX for VAR Modeling CCFTEST - Display CCF Function of Prewhitened data CHAR1 - Place a string is a character*1 array. CHARACTER - Place a string in a character*1 array. CHECKPOINT - Save workspace in portable file. CLEARALL - Clears all objects from workspace. CLEARDAT - Clears data from workspace. CLOSE - Close a logical unit. CLS - Clear screen. CLUSTER - K-Means and Hierarchical Cluster Models CMAXF1 - Constrained maximization of function using zxmwd. CMAXF2 - Constrained maximization of function using dbconf/g. CMAXF3 - Constrained maximization of function using db2pol. COMPRESS - Compress workspace. CONSTRAIN - Subset data based on range of values. CONTRACT - Contract a character array. COPY - Copy an object to another object COPYLOG - Copy file to log file. COPYOUT - Copy file to output file. COPYTIME - Copy time info from series 1 to series 2 COPYF - Copy a file from one unit to another. CSPECTRAL - Do cross spectral analysis. CSUB - Call Subroutine CSV - Read and Write a CSV file DATA_ACF - Calculate ACF and PACF Plots DATA2ACF - Calculate ACF and PACF Plots added argument DATAFREQ - Data Frequency DATAVIEW - View a Series Under Menu Control DELETECOL - Delete a column from a matrix or array. DELETEROW - Delete a row from a matrix or array. DES - Code / decode. DESCRIBE - Calculate Moment 1-4 and 6 of a series DF - Calculate Dickey-Fuller Unit Root Test. DISPLAYB - Displays a Buffer contents DIST_TAB - Distribution Table DMFGET - Load data froma DMF Save file DMFPUT - Place data in a DMF save file DMFMERGE - Merge two DMF save files DODOS - Execute a command string if under dos/windows. DO_SPEC - Display Periodogram and Spectrum DO2SPEC - Display Periodogram and Spectrum added argument DOUNIX - Execute a command string if under unix. DQDAG - Integrate a function using Gauss-Kronrod rules DQDNG - Integrate a smooth function using a nonadaptive rule. DQDAGI - Integrate a function over infinite/semi-infinite interval. DQDAGP - Integrate a function with singularity points given DQDAGS - Integrate a function with end point singularities DQAND - Multiple integration of a function DTWODQ - Two Dimensional Iterated Integral ESACF - Extended Sample Autocorrelation Function ECHOOFF - Turn off listing of execution. ECHOON - Turn on listing of execution. EPPRINT - Print to log and output file. EPRINT - Print to log file. ERASE - Erase file(s). EXPAND - Expand a character array FORMS - Build Control Forms FORPLOT - Forecast Plot FREE - Free a variable. FPLOT - Plot a Function FPRINT - Formatted print facility. GAMFIT - Generalized Additive Model Estimation GARCH - Calculate function for a ARCH/GARCH model. GARCHEST - Estimate ARCH/GARCH model. GET - Gets a variable from b34s. GETDMF - Gets a data from a b34s DFM file. GETKEY - Gets a key GETMATLAB - Gets data from matlab. GET_FILE - Gets a File name GET_NAME - Get Name of a Matrix Variable GETRATS - Reads RATS Portable file. GETSCA - Reads SCA FSAVE and MAD portable files. GMFAC - LU factorization of n by m matrix GMINV - Inverse of General Matrix using LAPACK GMSOLV - Solve Linear Equations system using LAPACK GRAPH - High Resolution graph. GRAPHP - Multi-Pass Graphics Programing Capability GRCHARSET - Set Character Set for Graphics. GRREPLAY - Graph replay and reformat command. GTEST - Tests output of a ARCH/GARCH Model GWRITE - Save Objects in GAUSS Format using one file GWRITE2 - Save objects in GAUSS format using two files HEADER - Turn on header HEXTOCH - Concert hex to a character representation. HINICH82 - Hinich 1982 Nonlinearity Test. HINICH96 - Hinich 1996 Nonlinearity Test. HPFILTER - Hodrick-Prescott Filter. ISEXTRACT - Place data in a structure. IALEN - Get actual length of a buffer of character data IBFCLOSE - Close a file that was used for Binary I/O IBFOPEN - Open a File for Binary I/O IBFREADC - Reads from a binary file into Character*1 array IBFREADR - Reads from a binary file into Real*8 array IBFSEEK - Position Binary read/write pointer IBFWRITER - Write noncharacter buffer on a binary file IBFWRITEC - Write character buffer on a binary file IB34S11 - Parse a token using B34S11 parser IFILESIZE - Determine number of bites in a file IFILLSTR - Fill a string with a character IGETICHAR - Obtain ichar info on a character buffer IGETCHARI - Get character from ichar value IJUSTSTR - Left/Right/center a string ILCOPY - Move bites from one location to another ILOCATESTR - Locate a substring in a string - 200 length max ILOWER - Lower case a string - 200 length max INEXTR8 - Convert next value in string to real*8 variable INEXTR4 - Convert next value in string to real*4 variable INEXTSTR - Extract next blank deliminated sub-string from a string INEXTI4 - Convert next value in a string to integer. INTTOSTR - Convert integer to string using format IRF - Impulse Response Functions of VAR Model IR8TOSTR - Convert real*8 value to string using format ISTRTOR8 - Convert string to real*8 ISTRTOINT - Convert string to integer IUPPER - Upper case a string - 200 length max I_DRNSES - Initializes the table used by shuffled generators. I_DRNGES - Get the table used in the shuffled generators. I_DRNUN - Uniform (0,1) Generator I_DRNNOR - Random Normal Distribution I_DRNBET - Random numbers from beta distribution I_DRNCHI - Random numbers from Chi-squared distribution I_DRNCHY - Random numbers from Cauchy distribution I_DRNEXP - Random numbers from standard exponential I_DRNEXT - Random numbers from mixture of two exponential distributions I_DRNGAM - Random numbers from standard gamma distribution I_DRNGCT - Random numbers from general continuous distribution I_DRNGDA - Random integers from discrete distribution alias approach I_DRNGDT - Random integers from discrete using table lookup I_DRNLNL - Random numbers from lognormal distribution I_DRNMVN - Random numbers from multivariate normal I_DRNNOA - Random normal numbers using acceptance/rejection I_DRNNOR - Random normal numbers using CDF method I_DRNSTA - Random numbers from stable distribution I_DRNTRI - Random numbers from triangular distribution I_DRNSPH - Random numbers on the unit circle I_DRNVMS - Random numbers from Von Mises distribution I_DRNWIB - Random numbers from Weibull distribution I_RNBIN - Random integers from binomial distribution I_RNGET - Gets seed used in IMSL Random Number generators. I_RNOPG - Gets the type of generator currently in use. I_RNOPT - Selects the type of uniform (0,1) generator. I_RNSET - Sets seed used in IMSL Random Number generators. I_RNGEO - Random integers from Geometric distribution I_RNHYP - Random integers from Hypergeometric distribution. I_RNMTN - Random numbers from multinomial distribution I_RNNBN - Negative binomial distribution I_RNPER - Random perturbation of integers I_RNSRI - Index of random sample without replacement KEENAN - Keenan Nonlinearity test KPSS - KPSS Stationarity Test KSWTEST - K Period Stock Watson Test KSWBOOTS - Critical values for output from KSWTEST KSWTESTM - Moving Period Stock Watson Test LAGMATRIX - Builds Lag Matrix. LAGTEST - 3-D Graph to display RSS for OLS Lags LAGTEST2 - 3-D Graph to display RSS for MARS Lags LAPACK - Sets Key LAPACK parameters LM - Engle Lagrange Multiplier ARCH test. LOAD - Load a Subroutine from a library. LOADDATA - Load Data from b34s into MATRIX command. LPMAX - Solve Linear Programming maximization problem. LPMIN - Solve Linear Programming minimization problem. LRE - McCullough Log Relative Error MAKEDATA - Place data in a b34s data loading structure. MAKEFAIR - Make Fair-Parke Data Loading File MAKEGLOBAL - Make a variable global (seen at all levels). MAKELOCAL - Make a variable seen at only local level. MAKEMATLAB - Place data in a file to be loaded into Matlab. MAKEMAD - Makes SCA *.MAD datafile from vectors MAKERATS - Make RATS portable file. MAKESCA - Make SCA FSV portable file. MANUAL - Place MATRIX command in manual mode. MARS - Multivariate Autoregressive Spline Models MARSPLINE - Updated MARS Command using Hastie-Tibshirani code MARS_VAR - Joint Estination of VAR Model using MARS Approach MAXF1 - Maximize a function using IMSL ZXMIN. MAXF2 - Maximize a function using IMSL DUMINF/DUMING. MAXF3 - Maximize a function using simplex method (DU2POL). MELD - Form all possible combinations of vectors. MENU - Put up user Menu for Input MESSAGE - Put up user message and allow a decision. MINIMAX - Estimate MINIMAX with MAXF2 MISSPLOT - Plot of a series with Missing Data MQSTAT - Multivariate Q Statistic MVNLTEST - Multivariate Third Order Hinich Test NAMES - List names in storage. NLEQ - Jointly solve a number of nonlinear equations. NLLSQ - Nonlinear Least Squares Estimation. NL2SOL - Alternative Nonlinear Least Squares Estimation. NLPMIN1 - Nonlinear Programming fin. diff. grad. DN2CONF. NLPMIN2 - Nonlinear Programming user supplied grad. DN2CONG. NLPMIN3 - Nonlinear Programming user supplied grad. DN0ONF. NLSTART - Generate starting values for NL routines. NOHEADER - Turn off header. OLSQ - Estimate OLS, MINIMAX and L1 models. OLSPLOT - Plot of Fitted and Actual Data & Res OPEN - Open a file and attach to a unit. OUTDOUBLE - Display a Real*8 value at a x, y on screen. OUTINTEGER - Display an Integer*4 value at a x, y on screen. OUTSTRING - Display a string value at a x, y point on screen. PCOPY - Copy an object from one pointer address to another PERMUTE - Reorder Square Matrix PISPLINE - Pi Spline Nonlinear Model Building PLOT - Line-Printer Graphics POLYFIT - Fit an nth degree polynomial POLYVAL - Evaluate an nth degree polynomial POLYMCONV - Convert storage of a polynomial matrix POLYMDISP - Display/Extract a polynomial matrix POLYMINV - Invert a Polynomial Matrix POLYMMULT - Multiply a Polynomial Matrix PPEXP - Two-dimensional exploratory projection pursuit PPEXP_P - Plot and save ppexp output PPREG - Projection Pursuit Regression PP - Calculate Phillips Peron Unit Root test PRINT - Print text and data objects. PRINTALL - Lists all variables in storage. PRINTOFF - Turn off Printing PRINTON - Turn on Printing (This is the default) PRINTVASV - Resets so that vectors/arrays print as vectors/arrays PRINTVASCMAT - Vectors/Arrays print as Column Matrix/Array PRINTVASRMAT - Vectors/Arrays print as Row Matrix/Array PROBIT - Estimate Probit (0-1) Model. PVALUE_1 - Present value of $1 recieved at end of n years PVALUE_2 - Present value of an Annuity of $1 PVALUE_3 - Present value of $1 recieved throughout year QPMIN - Quadratic Programming. QUANTILE - Calculate interquartile range. RANFOREST - Call the Random Forest Capability RCOVER - Recursive Covering DART/HYESS READ - Read data directly into MATRIX workspace from a file. REAL16INFO - Obtain Real16 info REAL16OFF - Turn off Real16 add REAL16ON - Turn on extended accuracy REAL32OFF - Turn off Real32 add REAL32ON - Turn on extended accuracy for real*16 REAL32_VPA - Turn on extended accuracy for real*16 using vpa RESET - Calculate Ramsey (1969) regression specification test. RESET77 - Thursby - Schmidt Regression Specification Test RESTORE - Load data back in MATRIX facility from external save file. RTEST - Test Residuals of Model RTEST2 - Test Residuals of Model - No RES and Y Plots REVERSE - Test a real*8 vector for reversibility in Freq. Domain REWIND - Rewind logical unit. ROTHMAN - Test a real*8 vector for reversibility in Time Domain RMATLAB - Runs Matlab RRPLOTS - Plots Recursive Residual Data RRPLOTS2 - Plots Recursive Residual Coef RUN - Terminates the matrix command being in "manual" mode. SAVE - Save current workspace in portable file format. SCHUR - Performs Schur decomposition SCREENCLOSE - Turn off Display Manager SCREENOPEN - Turn on Display Manager SCREENOUTOFF - Turn screen output off. SCREENOUTON - Turn screen output on. SET - Set all elements of an object to a value. SETCOL - Set column of an object to a value. SETLABEL - Set the label of an object. SETLEVEL - Set level. SETNDIMV - Sets an element in an n dimensional object. SETROW - Set row of an object to a value. SETTIME - Sets the time info in an existing series SETWINDOW - Set window to main(1), help(2) or error(3). SIGD - Set print digits. Default g16.8 SIMULATE - Dynamically Simulate OLS Model SMOOTH - Do exponential smoothing. SOLVEFREE - Set frequency of freeing temp variables. SORT - Sort a real vector. SPECTRAL - Spectral analysis of a vector or 1d array. STEPWISE - Stepwise OLS Regression STOP - Stop execution of a matrix command. SUBRENAME - Internally rename a subroutine. SUSPEND - Suspend loading and Execuiting a program SWARTEST - Stock-Watson VAR Test SYSTEM - Issue a system command. TABULATE - List vectors in a table. TB_TO_JULIAN - Converts tbase tstart freq to a real*8 julian date TESTARG - Lists what is passed to a subroutine or function. TIMER - Gets CPU time. TRIPLES - Calculate Triples Reversability Test TSAY - Calculate Tsay nonlinearity test. TSLINEUP - Line up Time Series Data TSD - Interface to TSD Data set UP_DATE_TB - Update the object header info UROOT_T - Battery of DF, PP, ERS and KPSS test VAREST - VAR Modeling VPASET - Set Variable Precision Math Options VOCAB - List built-in subroutine vocabulary. WRITE - Write an object to an external file. +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ Matrix Command Built-In Function Vocabulary +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ACF - Calculate autocorrelation function of a 1d object. AFAM - Change a matrix or vector to an array class object. ARGUMENT - Unpack character argument at run-time ARRAY - Define a 1d or 2d array. BETAPROB - Calculate a beta probability. BINDF - Evaluate Binomial Distribution Function BINPR - Evaluate Binomial Probability Function BOOTI - Calculate integers to be used with bootstrap. BOOTV - Bootstraps a vector with replacement. BOXCOX - Box-Cox Transformation of a series given lamda. BSNAK - Compute Not a Knot Sequence BSOPK - Compute optimal spline knot sequence BSINT - Compute 1-D spline interpolant given knots BSINT2 - Compute 2-D spline interpolant given knots BSINT3 - Compute 3-D spline interpolant given knots BSDER - Compute 1-D spline values/derivatives given knots BSDER2 - Compute 2-D spline values/derivatives given knots BSDER3 - Compute 3-D spline values/derivatives given knots BSITG - Compute 1-D spline integral given knots BSITG2 - Compute 2-D spline integral given knots BSITG3 - Compute 3-D spline integral given knots C1ARRAY - Create a Character*1 array C8ARRAY - Create a Character*8 array CATCOL - Concatenates an object by columns. CATROW - Concatenates an object by rows. CCF - Calculate the cross correlation function on two objects. CHAR - Convect an integer in range 0-127 to character. CHARDATE - Convert julian variable into character date dd\mm\yy. CHARDATEMY - Convert julian variable into character data mm\yyyy. CHARTIME - Converts julian variable into character date hh:mm:ss CHISQPROB - Calculate chi-square probability. CHTOR - Convert a character variable to a real variable. CHTOHEX - Convert a character to its hex representation. CFUNC - Call Function COMB - Combination of N objects taken M at a time. COMPLEX - Build a complex variable from two real*8 variables. CSPLINEFIT - Fit a 1 D Cubic Spline using alternative models CSPLINE - Calculate a cubic spline for 1 D data CSPLINEVAL - Calculate spline value given spline CSPLINEDER - Calculate spline derivative given spline value CSPLINEITG - Calculate integral of a cubic spline CUSUM - Cumulative sum. CUSUMSQ - Cumulative sum squared. CWEEK - Name of the day in character. DABS - Absolute value of a real*8 variable. DARCOS - Arc cosine of a real*8 variable. DARSIN - Arc sine of a real*8 variable. DATAN - Arc tan of a real*8 variable. DATAN2 - Arc tan of x / y. Signs inspected. DATENOW - Date now in form dd/mm/yy DBLE - Convert real*4 to real*8. DCONJ - Conjugate of complex argument. DCOS - Cosine of real*8 argument. DCOSH - Hyperbolic cosine of real*8 argument. DDOT - Inner product to two vectors. DERF - Error function of real*8/real*16 argument. DERFC - Inverse of error function. DERIVATIVE - Analytic derivative of a vector. DET - Determinate of a matrix. DEXP - Exponential of a real*8 argument. DFLOAT - Convert integer*4 to real*8. DGAMMA - Gamma function of real*8 argument. DIAG - Place diagonal of a matrix in an array. DIAGMAT - Create diagonal matrix. DIF - Difference a series. DINT - Extract integer part of real*8 number DNINT - Extract nearest integer part of real*8 number DIVIDE - Divide with an alternative return. DLGAMMA - Natural log of gamma function. DLOG - Natural log. DLOG10 - Base 10 log. DMAX - Largest element in an array. DMAX1 - Largest element between two arrays. DMIN - Smallest element in an array. DMIN1 - Smallest element between two arrays. DMOD - Remainder. DROPFIRST - Drops observations on top or array. DROPLAST - Drops observations on bottom of an array. DSIN - Calculates sine. DSINH - Hyperbolic sine. DSQRT - Square root of real*8 or complex*16 variable. DTAN - Tangent. DTANH - Hyperbolic tangent. EIGENVAL - Eigenvalue of matrix. Alias EIG. EPSILON - Positive value such that 1.+x ne 1. EVAL - Evaluate a character argument EXP - Exponential of real*8 or complex*16 variable. EXTRACT - Extract elements of a character*1 variable. FACT - Factorial FDAYHMS - Gets fraction of a day. FFT - Fast fourier transform. FIND - Finds location of a character string. FLOAT - Converts integer*4 to real*4. FPROB - Probability of F distribution. FREQ - Gets frequency of a time series. FRACDIF - Fractional Differencing FYEAR - Gets fraction of a year from julian date. GENARMA - Generate an ARMA series given parameters. GETDAY - Obtain day of year from julian series. GETHOUR - Obtains hour of the day from julian date. GETNDIMV - Obtain value from an n dimensional object. GETMINUTE - Obtains minute of the day from julian date. GETMONTH - Obtains month from julian date. GETQT - Obtains quarter of year from julian date. GETSECOND - Obtains second from julian date. GETYEAR - Obtains year. GOODCOL - Deletes all columns where there is missing data. GOODROW - Deletes all rows where there is missing data. GRID - Defines a real*8 array with a given increment. HUGE - Largest number of type HYPDF - Evaluate Hypergeometric Distribution Function HYPPR - Evaluate Hypergeometric Probability Function INTEGER8 - Load an Integer*8 object from a string I4TOI8 - Move an object from integer*4 to integer*8 I8TOI4 - Move an object from integer*8 to integer*4 ICHAR - Convect a character to integer in range 0-127. ICOLOR - Sets Color numbers. Used with Graphp. IDINT - Converts from real*8 to integer*4. IDNINT - Converts from real*8 to integer*4 with rounding. INFOGRAPH - Obtain Interacter Graphics INFO IMAG - Copy imaginary part of complex*16 number into real*8. IAMAX - Largest abs element in 1 or 2D object IAMIN - Smallest abs element in 1 or 2D object IMAX - Largest element in 1 or 2D object IMIN - Smallest element in 1 or 2D object INDEX - Define integer index vector, address n dimensional object. INLINE - Inline creation of a program INT - Copy real*4 to integer*4. INTEGERS - Generate an integer vector with given interval. INV - Inverse of a real*8 or complex*16 matrix. INVBETA - Inverse beta distribution. INVCHISQ - Inverse Chi-square distribution. INVFDIS - Inverse F distribution. INVTDIS - Inverse t distribution. IQINT - Converts from real*16 to integer*4. IQNINT - Converts from real*16 to integer*4 with rounding. ISMISSING - Sets to 1.0 if variable is missing IWEEK - Sets 1. for monday etc. JULDAYDMY - Given day, month, year gets julian value. JULDAYQY - Given quarter and year gets julian value. JULDAYY - Given year gets julian value. KEEPFIRST - Given k, keeps first k observations. KEEPLAST - Given k, keeps last k observations. KIND - Returns kind of an object in integer. KINDAS - Sets kind of second argument to kind first arg. KLASS - Returns klass of an object in integer. KPROD - Kronecker Product of two matrices. LABEL - Returns label of a variable. LAG - Lags variable. Missing values propagated. LEVEL - Returns current level. LOWERT - Lower triangle of matrix. MCOV - Consistent Covariance Matrix MAKEJUL - Make a Julian date from a time series MASKADD - Add if mask is set. MASKSUB - Subtract if mask is set. MATRIX - Define a matrix. MEAN - Average of a 1d object. MEDIAN - Median of a real*8 object MFAM - Set 1d or 2d array to vector or matrix. MISSING - Returns missing value. MLSUM - Sums log of elements of a 1d object. MOVELEFT - Moves elements of character variable left. MOVERIGHT - Move elements of character variable right. NAMELIST - Creates a namelist. NEAREST Nearest distinct number of a given type NCCHISQ - Non central chi-square probability. NOCOLS - Gets number of columns of an object. NOELS - Gets number of elements in an object. NORMDEN - Normal density. NORMDIST - 1-norm, 2-norm and i-norm distance. NOROWS - Gets number of rows of an object. NOTFIND - Location where a character is not found. OBJECT - Put together character objects. PDFAC - Cholesky factorization of PD matrix. PDFACDD - Downdate Cholesky factorization. PDFACUD - Update Cholesky factorization. PDINV - Inverse of a PD matrix. PDSOLV - Solution of a PD matrix given right hand side. PI - Pi value. PINV - Generalized inverse. PLACE - Places characters inside a character array. POIDF - Evaluate Poisson Distribution Function POIPR - Evaluate Poisson Probability Function POINTER - Machine address of a variable. POLYDV - Division of polynomials. POLYMULT - Multiply two polynomials POLYROOT - Solution of a polynomial. PROBIT - Inverse normal distribution. PROBNORM - Probability of normal distribution. PROBNORM2 - Bivariate probability of Nornal distribution. PROD - Product of elements of a vector. Q1 - Q1 of a real*8 object Q3 - Q3 of a real*8 object QCOMPLEX - Build complex*32 variable from real*16 inputs. QINT - Extract integer part of real*16 number QNINT - Extract nearest integer part of real*16 number QREAL - Obtain real*16 part of a complex*326 number. QR - Obtain Cholesky R via QR method using LAPACK. QRFAC - Obtain Cholesky R via QR method. QRSOLVE - Solve OLS using QR. RANKER - Index array that ranks a vector. RCOND - 1 / Condition of a Matrix REAL - Obtain real*8 part of a complex*16 number. R8TOR16 - Convert Real*8 to Real*16 R16TOR8 - Convert Real*16 to Real*8 REAL16 - Input a Real*16 Variable REC - Rectangular random number. RECODE - Recode a real*8 or character*8 variable RN - Normally distributed random number. ROLLDOWN - Moves rows of a 2d object down. ROLLLEFT - Moves cols of a 2d object left. ROLLRIGHT - Moves cols of a 2d object right. ROLLUP - Moves rows of a 2d object up. RTOCH - Copies a real*8 variable into character*8. SEIGENVAL - Eigenvalues of a symmetric matrix. Alias SEIG. SEXTRACT - Takes data out of a field. SFAM - Creates a scalar object. SNGL - Converts real*8 to real*4. SPACING - Absolute spacing near a given number SPECTRUM - Returns spectrum of a 1d object. SUBSET - Subset 1d, 2d array, vector or matrix under a mask. SUBMATRIX - Define a Submatrix SUM - Sum of elements. SUMCOLS - Sum of columns of an object. SUMROWS - Sum of rows of an object. SUMSQ - Sum of squared elements of an object. SVD - Singular value decomposition of an object. TIMEBASE - Obtains time base of an object. TIMENOW - Time now in form hh:mm:ss TIMESTART - Obtains time start of an object. TINY - Smallest number of type TDEN - t distribution density. TO_RMATRIX - Convert Object to Row-Matrix TO_CMATRIX - Convert Object to Col-Matrix TO_RARRAY - Convert Object to Row-Array TO_CARRAY - Convert Object to Col-Matrix TO_VECTOR - Convert Object to Vector TO_ARRAY - Convert Object to Array TPROB - t distribution probability. TRACE - Trace of a matrix. TRANSPOSE - Transpose of a matrix. UPPERT - Upper Triangle of matrix. VARIANCE - Variance of an object. VECTOR - Create a vector. VFAM - Convert a 1d array to a vector. VOCAB - List built in functions. VPA - Variable Precision Math calculation ZDOTC - Conjugate product of two complex*16 objects. ZDOTU - Product of two complex*16 objects. ZEROL - Zero lower triangle. ZEROU - Zero upper triangle. ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ Matrix Programming Language key words ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ CALL - Call a subroutine CONTINUE - go to statement DO - Starts a do loop DOWHILE - Starts a dowhile loop NEXT i - End of a do loop ENDDO - End of a do loop ENDDOWHILE - End of a dowhile loop END - End of a program, function or Subroutine. EXITDO - Exit a DO loop EXITIF - Exit an IF statement FOR - Start a do loop, FORMULA - Define a recursive formula. GO TO - Transfer statement FUNCTION - Beginning of a function. IF( ) - Beginning of an IF structure ENDIF - End of an IF( )THEN structure PROGRAM - Beginning of a program, RETURN - Next to last statement before end. RETURN( ) - Returns the result of a function. SOLVE - Solve a recursive system. SUBROUTINE - Beginning of subroutine. WHERE( ) - Starts a where structure. ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ matrix2.mac capability as of 1 May 2007 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ Mac file name c:\b34slm\lib\matrix2.mac Name Heading ACF_PLOT Plot ACF ALTSE Copies lower triangle to upper and gets SE AUTOCOV Autocovaiance BJ_IDEN List / Plot ACF and PACF B_G_TEST Breusch- Godfrey (1978) Test on Residuals B_G_ALT Breusch- Godfrey (1978) Test using dropping BLUS BLUS master BPF Baxter - King Symmetric MA Filter BPFM Baxter - King Symmetric MA Filter - Missing Data at BUILDLAG Builds NEWY and NEWX for VAR Modeling CCFTEST Tests Cross Correlations CFREQ Determine Frequency Distribution COINT2 Cointegreation of Two Series COINT2LM Cointegreation of Two Series - Gives L1 and Minimax COINT2M Moving Cointegration of Two Series - Simple Version COINT2M2 Moving Cointegration Two Series OLS, L1 Minimax COINT2ME Moving Cointegration of Two Series - Extended Versi COINT3 Cointegration of Three Series COINT3ME Moving Cointegration of Three Series DATAVIEW Interactively View Data DATA_ACF Plot ACF and PACF of a series DATA2ACF Plot ACF and PACF of a series with file name DF_GLS Elliot-Rothenberg-Stock DF_GLS Test DIST_TAB Distribution Table DO_SPEC View Spectrum DO2SPEC View Spectrum with a file name DUD Derivative Free Nonlinear Estimation DUD2 Derivative Free Nonlinear Estimation FDIFINFO Fractional Differencing FILTER High Low Pass Filter using Real FFT FILTERC High Low Pass Filter using Complex FFT FORPLOT Forecast Plot GAMFORE Forecast a GAM Model GAMPLOT Display GAMFIT Results GARCH2P Two Pass GARCH Using ARMA Command GARCH2PF Two Pass GARCH Using ARMA Command with Forecasts GET_NAME Get Name of a Matrix Variable GTEST Tests ARCH / GARCH Models GWRITE Save Matrix Object in GAUSS Format using one file GWRITE2 Pass Large Datasets to GAUSS in 2 files HP_2 Hodrick - Prescott Moving Filtering HP_BP_1 Baxter - King & Hodrick - Prescott Filtering HP_BP_2 Baxter - King & Hodrick - Prescott Moving Filtering IRF Impulse Response Functions for VAR Model KPSS KPSS Stationarity Test KSWTEST K Period Stock Watson Test KSWBOOTS Critical values for output from KSWTEST KSWTESTM Moving Period Stock Watson Test LAGTEST 3-D Graph to display RSS for Various OLS Lags LAGTEST2 3-D Graph to display RSS for Various MARS Lags LMTEST Engle LM test for a Range of lags MARQ Estimation of a Nonlinear Model using Derivatives MARSPLOT Plot MARS Curves & Surfaces MCLEODLI McLeod-Li (1983) Linearity test (y,ip,maxacf) MINIMAX Minimax Estimation using MAXF2 MISSPLOT Plots Data With Missing Values MOVEAVE Moving average of a vector MOVEBJ Moving Arima Forecast using AUTOBJ MOVECORR Moving Correlation of two vectors MOVEH82 Moving Hinich 82 test MOVEH96 Moving Hinich 96 test MOVEOLS Moving OLS with LAGS MOVEVAR Moving Variance NLVARCOV NLLSQ Variance Covariance OLSPLOT Plot of Fitted and Actual Data & Res PAD Pad a 1D Real*8 Series on both ends PERMUTE Reorder a Square Matrix POLYFIT Fit an nth Order Regression POLYVAL Evaluate an nth Order Polynominam Regression PVALUE_1 Present value of $1 recieved at end of n years PVALUE_2 Present Value of an Annuity of $1 PVALUE_3 Present value of $1 recieved throughout year QUANTREG Quantile Regression RESET69 Ramsey(1969) Regression Specification Test RESET77 RESET77 Nonlinearitry Test RMATLAB Runs Matlab commands RRPLOTS Plots Recursive Residual Data RRPLOTS2 Plots Recursive Coef and t scores RTEST Tests Residuals of TS Models RTEST2 Tests Residuals of TS Models - No Res and y Plots SUBSET Subset an array under mask control SWARTEST Stock-Watson AR Test TESTFUN Test Function TESTPGM Test Program TESTSUB Test Subroutine UROOT_T Battery of DF, PP, ERS and KPSS test VAREST VAR Modeling ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ staging2.mac capability as of 1 May 2007 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ Mac file name c:\b34slm\lib\staging2.mac Name Heading ACE_OLS Tests ACE Model using OLS ACE_PLOT Plot ACE Response Curves BOOTOLS Bootstrap OLS Model BOOTL1 Bootstrap L1 Model BOOTMM Bootstrap MM Model COR Correlation of a matrix. COV Calculated Covariance of a matrix COND Condition of a Matrix DFVALUES Estimate of DF/PP Critical Values DISPMARS Displays Mars Model F_MARSLG Adjusts yhat and forecasts from MARSPLINE G4TEST Pena-Slate JASA March (2006) Joint OLS Test GAMPLOT Display GAMFIT Results GARCHF Forecast ARCH / GARCH Series GARCH2PA Automatic GARCH Two Pass Modeling GENGARCH Simulate GARCH and ARCH Models GETCHAR Unpack a string into Character Array GETDATA Simple Real*8 data Read routine GETDATA2 Simple Real*16 data Read routine GETR16 Unpack a string into Real*16 Array GETVPA Unpack a string into VPA Array GETR8 Unpack a string into Real*8 Array GLESJER Glesjer (1969) Heteroscedasticity Test GLS Cochrane & Orcutt GLS Estimation of AR(k) Model GLSDATA Transform x to x* given GLS rho vector GLS_ML GLS_ML Maximum Likelihood estimation of AR(1) Model GLS_2P Two Pass Generalized Least Squares for Order K GPH Geweke & Porter-Hudak (1983) Fractional Diff. Test G_QUANDT Goldfeld-Quandt Heteroscedasticity Test HET_TEST Heteroscedasticity Testing IACF Inverse ACF LASSO LASSO Shrinkage Approach MARSINFO Estimate of MARSPLINE Model GCV Info MCOVF Function for Covariance. See built in mcov( ) NORM NORM(X,i) of a matrix or vector NTOKIN Counts number of tokins in a string NW_SE Newey West SE Correction Calculation OLS_CON Constrained OLS Model P_L_EST Probit & Logit Estimation PANEL_LIB Panel Subroutine Library QR_SMALL Break Appart QR Factorization RIDGE Ridge Regression RNW_SE Run Newey West SE Correction for OLS Model SC_TEST Serial Correlation Tests SVD_OLS OLS Using SVD Approach SVD2_OLS OLS Using SVD Approach and LAPACK TLOGIT Tests Logit/Probit Model TS_TESTS Residual Specification tests TSALIGN Align Time Series Data UNDIF Undifference a series WALD Wald Test of restrictions +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ Commands listed by Function Character & String Subroutines CHTOHEX - Convert a character to its hex representation. CONTRACT - Contract a character array DISPLAYB - Displays a Buffer contents EXPAND - Expand a character array EVAL - Evaluate a Character Argument HEXTOCH - Concert hex to a character representation. IALEN - Get actual length of a buffer of character data IBFCLOSE - Close a file that was used for Binary I/O IBFOPEN - Open a File for Binary I/O IBFREADC - Reads from a binary file into Character*1 array IBFREADR - Reads from a binary file into Real*8 array IBFSEEK - Position Binary read/write pointer IBFWRITER - Write noncharacter buffer on a binary file IBFWRITEC - Write character buffer on a binary file IB34S11 - Parse a token using B34S11 parser IFILESIZE - Determine number of bites in a file IFILLSTR - Fill a string with a character IGETICHAR - Obtain ichar info on a character buffer IGETCHARI - Get character from ichar value IJUSTSTR - Left/Right/center a string ILCOPY - Move bites from one location to another ILOCATESTR - Locate a substring in a string - 200 length max ILOWER - Lower case a string - 200 length max INEXTR8 - Convert next value in string to real*8 variable INEXTR4 - Convert next value in string to real*4 variable INEXTSTR - Extract next blank deliminated sub-string from a string INEXTI4 - Convert next value in a string to integer. INTTOSTR - Convert integer to string using format IR8TOSTR - Convert real*8 value to string using format ISTRTOR8 - Convert string to real*8 ISTRTOINT - Convert string to integer IUPPER - Upper case a string - 200 length max Character Functions CHAR - Convect an integer in range 0-127 to character. CHTOR - Convert a character variable to a real variable. DATENOW - Date now in form dd/mm/yy EXTRACT - Extract elements of a character*1 variable. FIND - Finds location of a character string. ICHAR - Convect a character to integer in range 0-127. INLINE - Inline creation of a program MOVELEFT - Moves elements of character variable left. MOVERIGHT - Move elements of character variable right. NAMELIST - Creates a namelist. NOTFIND - Location where a character is not found. OBJECT - Put together character objects. PLACE - Places characters inside a character array. RTOCH - Copies a real*8 variable into character*8. TIMENOW - Time now in form hh:mm:ss Data Building, loading and Display Subroutines ACF_PLOT - Simple ACF Plot ADDCOL - Add a column to a 2d array or matrix. ADDROW - Add a row to a 2d array or matrix. AGGDATA - Aggregate Data under control of an ID Vector. ALIGN Align Series with Missing Data BJ_IDEN - ACF and PACF Listing and or Plotting BUILDLAG - Builds NEWY and NEWX for VAR Modeling CCFTEST - Display CCF Function of Prewhitened data CHAR1 - Place a string in a character*1 array. CHARACTER - Place a string in a character*1 array. CHECKPOINT - Save workspace in portable file. CLEARALL - Clears all objects from workspace. CLEARDAT - Clears data from workspace. CSUB - Call Subroutine CONSTRAIN - Subset data based on range of values. CSV - Read and Write a CSV file DATA_ACF - Calculate ACF and PACF Plots DATA2ACF - Calculate ACF and PACF Plots added argument DATAVIEW - View a Series Under Menu Control DATAFREQ - Data Frequency DELETECOL - Delete a column from a matrix or array. DELETEROW - Delete a row from a matrix or array. DES - Code / decode. DESCRIBE - Calculate Moment 1-4 and 6 of a series DIVIDE - Divide with an alternative return. DO_SPEC - Display Periodogram and Spectrum DO2SPEC - Display Periodogram and Spectrum added argument ESACF - Extended Sample Autocorrelation Function GET_FILE - Get a File Name GET_NAME - Get Name of a Matrix Variable GTEST - Tests output of a ARCH/GARCH Model GWRITE - Save Objects in GAUSS Format using one file GWRITE2 - Save Objects in GAUSS Format using two files FREE - Free a variable. LAGMATRIX - Builds Lag matrix. MELD - Form all possible combinations of vectors. QUANTILE - Calculate interquartile range. RTEST - Test Residuals of Model RTEST2 - Test Residuals of OLS Model - No RES and Y Plots SET - Set all elements of an object to a value. SETCOL - Set column of an object to a value. SETLABEL - Set the label of an object. SETNDIMV - Sets an element in an n dimensional object. SETROW - Set row of an object to a value. SETTIME - Sets the time info in an existing series SORT - Sort a real vector. SUBRENAME - Internally rename a subroutine. Data Building Functions AFAM - Change a matrix or vector to an array class object. ARRAY - Define a 1d or 2d array. C1ARRAY - Create a Character*1 array C8ARRAY - Create a Character*8 array CATCOL - Concatenates an object by columns. CATROW - Concatenates an object by rows. CFUNC - Call Function COMB - Combination of N objects taken M at a time. COMPLEX - Build a complex variable from two real*8 variables. CUSUM - Cumulative sum. CUSUMSQ - Cumulative sum squared. DABS - Absolute value of a real*8 variable. DARCOS - Arc cosine of a real*8 variable. DARSIN - Arc sine of a real*8 variable. DATAN - Arc tan of a real*8 variable. DATAN2 - Arc tan of x / y. Signs inspected. DBLE - Convert real*4 to real*8. DCONJ - Conjugate of complex argument. DCOS - Cosine of real*8 argument. DCOSH - Hyperbolic cosine of real*8 argument. DDOT - Inner product to two vectors. DEXP - Exponential of a real*8 argument. DFLOAT - Convert integer*4 to real*8. DGAMMA - Gamma function of real*8 argument. DINT - Extract integer part of real*8 number DNINT - Extract integer part of real*8 number DLGAMMA - Natural log of gamma function. DLOG - Natural log. DLOG10 - Base 10 log. DMAX - Largest element in an array. DMAX1 - Largest element between two arrays. DMIN - Smallest element in an array. DMIN1 - Smallest element between two arrays. DMOD - Remainder. DROPFIRST - Drops observations on top or array. DROPLAST - Drops observations on bottom of an array. DSIN - Calculates sine. DSINH - Hyperbolic sine. DSQRT - Square root of real*8 or complex*16 variable. DTAN - Tangent. DTANH - Hyperbolic tangent. EPSILON - Positive value such that 1.+x ne 1. EXP - Exponential of real*8 or complex*16 variable. FLOAT - Converts integer*4 to real*4. FREQ - Gets frequency of a time series. FYEAR - Gets fraction of a year from julian date. GETIDIM - Obtain value from an n dimensional object. GMFAC - LU factorization of n by m matrix GMINV - Inverse of General Matrix using LAPACK GMSOLV - Solve Linear Equations system using LAPACK GOODCOL - Deletes all columns where there is missing data. GOODROW - Deletes all rows where there is missing data. GRID - Defines a real*8 array with a given increment. HUGE - Largest number of type IAMAX - Largest abs element in 1 or 2D object IAMIN - Smallest abs element in 1 or 2D object IMAX - Largest element in 1 or 2D object IMIN - Smallest element in 1 or 2D object INTEGER8 - Load an Integer*8 object from a string I4TOI8 - Move an object from integer*4 to integer*8 I8TOI4 - Move an object from integer*8 to integer*4 ICOLOR - Sets Color numbers. Used with Graphp. IDINT - Converts from real*8 to integer*4. IDNINT - Converts from real*8 to integer*4 with rounding. IMAG - Copy imaginary part of complex*16 number into real*8. INDEX - Define integer index vector. INFOGRAPH - Obtain Interacter Graphics INFO INT - Copy real*4 to integer*4. INTEGERS - Generate an integer vector with given interval. IQINT - Converts from real*16 to integer*4. IQNINT - Converts from real*16 to integer*4 with rounding. ISMISSING - Sets to 1.0 if variable is missing KEEPFIRST - Given k, keeps first k observations. KEEPLAST - Given k, keeps last k observations. KIND - Returns kind of an object in integer. KINDAS - Sets kind of second argument to kind first arg. KLASS - Returns klass of an object in integer. LABEL - Returns label of a variable. LAG - Lags variable. Missing values propagated. LOWERT - Lower triangle of matrix. MCOV - Consistent Covariance Matrix MASKADD - Add if mask is set. MASKSUB - Subtract if mask is set. MATRIX - Define a matrix. MFAM - Set 1d or 2d array to vector or matrix. MISSING - Returns missing value. NEAREST Nearest distinct number of a given type NOCOLS - Gets number of columns of an object. NOELS - Gets number of elements in an object. NOROWS - Gets number of rows of an object. NORMDIST - 1-norm, 2-norm and i-norm distance. PI - Pi value. QCOMPLEX - Build complex*32 variable from real*16 inputs. QINT - Extract integer part of real*16 number QNINT - Extract nearest integer part of real*16 number QREAL - Obtain real*16 part of a complex*32 number. RANKER - Index array that ranks a vector. REAL - Obtain real*8 part of a complex*16 number. R8TOR16 - Convert Real*8 to Real*16 R16TOR8 - Convert Real*16 to Real*8 REAL16 - Input a Real*16 Variable REC - Rectangular random number. RECODE - Recode a real*8 or chartacter*8 variable RN - Normally distributed random number. SEXTRACT - Takes data out of a field. SFAM - Creates a scalar object. SNGL - Converts real*8 to real*4. SPACING Absolute spacing near a given number SUBMATRIX - Define a Submatrix SUBSET - Subset 1d, 2d array, vector or matrix under a mask. SUM - Sum of elements. SUMCOLS - Sum of columns of an object. SUMROWS - Sum of rows of an object. SUMSQ - Sum of squared elements of an object. TIMEBASE - Obtains time base of an object. TIMESTART - Obtains time start of an object. TINY Smallest number of type TO_RMATRIX - Convert Object to Row-Matrix TO_CMATRIX - Convert Object to Col-Matrix TO_RARRAY - Convert Object to Row-Array TO_CARRAY - Convert Object to Col-Matrix TO_VECTOR - Convert Object to Vector TO_ARRAY - Convert Object to Array UPPERT - Upper Triangle of matrix. VECTOR - Create a vector. VFAM - Convert a 1d array to a vector. ZDOTC - Conjugate product of two complex*16 objects. ZDOTU - Product of two complex*16 objects. ZEROL - Zero lower triangle. ZEROU - Zero upper triangle. Data Filtering Subroutines CSPECTRAL - Do cross spectral analysis. SPECTRAL - Spectral analysis of a vector or 1d array. LOESS_MV - Multivariate Loess Smoothing LOESS_SA - Seasonal Adjustment using LOESS Method. LOESS_SPS - Scatter Plot Smoothing using LOESS Method. FFT - Fast fourier transform. PPEXP - Two-dimensional exploratory projection pursuit PPEXP_P - Plot and save ppexp output PPREG - Projection Pursuit Regression RANFOREST - Call the Random Forest Capability RCOVER - Recursive Covering DART/HYESS SPECTRUM - Returns spectrum of a 1d object. VAREST - VAR Modeling WAVELET - Calculate Wavelets Date and Time Functions CHARDATE - Convert julian variable into character date dd\mm\yy. CHARDATEMY - Convert julian variable into character data mm\yyyy. CHARTIME - Converts julian variable into character date hh:mm:ss COPYTIME - Copy time info from series 1 to series 2 CWEEK - Name of the day in character. FDAYHMS - Gets fraction of a day. GETDAY - Obtain day of year from julian series. GETHOUR - Obtains hour of the day from julian date. GETMINUTE - Obtains minute of the day from julian date. GETMONTH - Obtains month from julian date. GETQT - Obtains quarter of year from julian date. GETSECOND - Obtains second from julian date. GETYEAR - Obtains year. IWEEK - Sets 1. for monday etc. JULDAYDMY - Given day, month, year gets julian value. JULDAYQY - Given quarter and year gets julian value. JULDAYY - Given year gets julian value. MAKEJUL - Make a Julian date from a time series Estimation and residual testing Subroutines ARMA - ARMA estimation using ML and MOM. AUTOBJ - Automatic Estimation of Box-Jenkins Model BDS - BDS Nonlinearity test. BESTREG - Best OLS REGRESSION BPFILTER - Baxter-King Filter. DF - Calculate Dickey-Fuller Unit Root Test. GAMFIT - Generalized Additive Model Estimation GAMFORE - Forecast a GAM Model GAMPLOT - Display GAMFIT Results GARCHEST - Estimate ARCH/GARCH model. IRF - Impulse Response Functions of VAR Model KEENAN - Keenan Nonlinearity test KSWTEST - K Period Stock Watson Test KSWBOOTS - Critical values for output from KSWTEST KSWTESTM - Moving Period Stock Watson Test LAGTEST - 3-D Graph to display RSS for OLS Lags LAGTEST2 - 3-D Graph to display Rss for MARS Lags LM - Engle Lagrange Multiplier ARCH test. MINIMAX - Estimate MINIMAX with MAXF2 OLSQ - Estimate OLS, MINIMAX and L1 models. OLSPLOT - Plot of Fitted and Actual Data & Res POLYFIT - Fit an nth degree polynomial POLYVAL - Evaluate an nth degree polynomial PP - Calculate Phillips Peron Unit Root test PROBIT - Estimate Probit (0-1) Model. QUANTREG - Quantile Regression Program RESET - Calculate Ramsey(1969) regression specification test. RESET77 - Thursby - Schmidt Regression Specification Test REVERSE - Test a real*8 vector for reversibility in Freq. Domain ROTHMAN - Test a real*8 vector for reversibility in Time Domain RRPLOTS - Plots Recursive Residual Data RRPLOTS2 - Plots Recursive Residual Coef SIMULATE - Dynamically Simulate OLS Model SMOOTH - Do exponential smoothing. SWARTEST - Stock-Watson VAR Test STEPWISE - Stepwise OLS Regression TRIPLES - Calculate Triples Reversability Test TSAY - Calculate Tsay nonlinearity test. TSD - Interface to TSD Data set Estimation and residual testing and Spline Functions BOOTI - Calculate integers to be used with bootstrap. BOOTV - Bootstraps a vector with replacement. BOXCOX - Box-Cox Transformation of a series given lamda. MLSUM - Sums log of elements of a 1d object. Spline Related Commands ACEFIT - Alternating Conditional Expectation Model Estimation BSNAK - Compute Not a Knot Sequence BSOPK - Compute optimal spline knot sequence BSINT - Compute 1-D spline interpolant given knots BSINT2 - Compute 2-D spline interpolant given knots BSINT3 - Compute 3-D spline interpolant given knots BSDER - Compute 1-D spline values/derivatives given knots BSDER2 - Compute 2-D spline values/derivatives given knots BSDER3 - Compute 3-D spline values/derivatives given knots BSITG - Compute 1-D spline integral given knots BSITG2 - Compute 2-D spline integral given knots BSITG3 - Compute 3-D spline integral given knots CSPLINEFIT - Fit a 1 D Cubic Spline using alternative models CSPLINE - Calculate a cubic spline for 1 D data CSPLINEVAL - Calculate spline value given spline CSPLINEDER - Calculate spline derivative given spline value CSPLINEITG - Calculate integral of a cubic spline GAMFIT - Generalized Additive Model Estimation GAMFORE - Forecast a GAM Model GAMPLOT - Display GAMFIT Results MARS - Multivariate Autoregressive Spline Models MARSPLINE - Updated MARS Command using Hastie-Tibshirani code MARS_VAR - Joint Estination of VAR Model using MARS Approach PISPLINE - Pi Spline Nonlinear Model Building Integration Subroutines DQDAG - Integrate a function using Gauss-Kronrod rules DQDNG - Integrate a smooth function using a nonadaptive rule. DQDAGI - Integrate a function over infinite/semi-infinite interval. DQDAGP - Integrate a function with singularity points given DQDAGS - Integrate a function with end point singularities DQAND - Multiple integration of a function DTWODQ - Two Dimensional Iterated Integral IO routines Subroutines BACKSPACE - Backspace a unit CLOSE - Close a logical unit. COPYLOG - Copy file to log file. COPYOUT - Copy file to output file. COPYF - Copy a file from one unit to another. DMFGET - Load data froma DMF Save file DMFPUT - Place data in a DMF save file DMFMERGE - Merge two DMF save files ECHOOFF - Turn off listing of execution. ECHOON - Turn on listing of execution. EPPRINT - Print to log and output file. EPRINT - Print to log file. ERASE - Erase file(s). FORMS - Build Control Forms FPRINT - Formatted print facility. GET - Gets a variable from b34s. GETDMF - Gets a data from a b34s DFM file. GETKEY - Gets a key GETMATLAB - Gets data from matlab. GETRATS - Reads RATS Portable file. GETSCA - Reads SCA FSAVE and portable portable files HEADER - Turn on header ISEXTRACT - Place data in a structure. LOADDATA - Load Data from b34s into MATRIX command. MAKEDATA - Place data in a b34s data loading structure. MAKEDMF - Place Data in a b34s DMF file MAKEFAIR - Make Fair-Parke Data Loading File MAKEMAD - Makes SCA *.MAD datafile from vectors MAKEMATLAB - Place data in a file to be loaded into Matlab. MAKERATS - Make RATS portable file. MAKESCA - Make SCA FSV portable file. MENU - Put up user Menu for input MESSAGE - Put up user message and allow a decision. NOHEADER - Turn off header OPEN - Open a file and attach to a unit. PRINT - Print text and data objects. PRINTALL - Lists all variables in storage. PRINTOFF - Turn off Printing PRINTON - Turn on Printing (This is the default) READ - Read data directly into MATRIX workspace from a file. RESTORE - Load data back in MATRIX facility from external save file. REWIND - Rewind logical unit. RMATLAB - Runs Matlab SAVE - Save current workspace in portable file format. TABULATE - List vectors in a table. WRITE - Write an object to an external file. Matrix Functions DERIVATIVE - Analytic derivative of a vector. DET - Determinate of a matrix. DIAG - Place diagonal of a matrix in an array. DIAGMAT - Create diagonal matrix. EIGENVAL - Eigenvalue of matrix. Alias EIG. INV - Inverse of a real*8 or complex*16 matrix. KPROD - Kronecker Product of two matrices. PDFAC - Cholesky factorization of PD matrix. PDFACDD - Downdate Cholesky factorization. PDFACUD - Update Cholesky factorization. PDINV - Inverse of a PD matrix. PDSOLV - Solution of a PD matrix given right hand side. PERMUTE - Reorder Square Matrix PINV - Generalized Inverse. PROD - Product of elements of a vector. QR - Obtain Cholesky R via QR method using LAPACK. QRFAC - Obtain Cholesky R via QR method. QRSOLVE - Solve OLS using QR. RCOND - 1 / Condition of a Matrix. ROLLDOWN - Moves rows of a 2d object down. ROLLLEFT - Moves cols of a 2d object left. ROLLRIGHT - Moves cols of a 2d object right. ROLLUP - Moves rows of a 2d object up. SCHUR - Performs Schur decomposition SEIGENVAL - Eigenvalues of a symmetric matrix. Alias SEIG. SVD - Singular value decomposition of an object. TRACE - Trace of a matrix. TRANSPOSE - Transpose of a matrix. ZDOTC - Conjugate product of two complex*16 objects. ZDOTU - Product of two complex*16 objects. Miscellaneous Subroutines BREAK - Set User Program Break Point. COMPRESS - Compress workspace. COPY - Copy an object to another object LAPACK - Sets Key LAPACK parameters LOAD - Load a Subroutine from a library. MAKEGLOBAL - Make a variable global (seen at all levels). MAKELOCAL - Make a variable seen at only local level. MANUAL - Place MATRIX command in manual mode. NAMES - List names in storage. LRE - McCullough Log Relative Error PCOPY - Copy an object from one pointer address to another REAL16INFO - Obtain Real16 info REAL16OFF - Turn off Real16 add REAL16ON - Turn on extended accuracy REAL32OFF - Turn off Real32 add REAL16ON - Turn on Real*32 extended accuracy RUN - Terminates the matrix command being in "manual" mode. SETLEVEL - Set level. SETWINDOW - Set window to main(1), help(2) or error(3). SIGD - Set print digits. Default g16.8 STOP - Stop execution of a matrix command. TESTARG - Lists what is passed to a subroutine or function. TIMER - Gets CPU time. VOCAB - List built-in subroutine vocabulary. Miscellaneous Functions ARGUMENT - Unpack character argument at run-time FACT - Factorial LEVEL - Returns current level. POINTER - Machine address of a variable. POLYROOT - Solution of a polynomial. POLYDV - Division of polynomials. POLYMULT - Multiply two polynomials VOCAB - List built in functions. Nonlinear Estimation Capability Subroutines CMAXF1 - Constrained maximization of function using zxmwd. CMAXF2 - Constrained maximization of function using dbconf/g. CMAXF3 - Constrained maximization of function using db2pol. BGARCH Calculate function for a BGARCH model. GARCH - Calculate function for a ARCH/GARCH model. GARCHEST - Estimate ARCH/GARCH model. LPMAX - Solve Linear Programming maximization problem. LPMIN - Solve Linear Programming minimization problem. MAXF1 - Maximize a function using IMSL ZXMIN. MAXF2 - Maximize a function using IMSL DUMINF/DUMING. MAXF3 - Maximize a function using simplex method (DU2POL). NLEQ - Jointly solve a number of nonlinear equations. NLLSQ - Nonlinear Least Squares Estimation. NL2SOL - Alternative Nonlinear Least Squares Estimation. NLPMIN1 - Nonlinear Programming fin. diff. grad. DN2CONF. NLPMIN2 - Nonlinear Programming user supplied grad. DN2CONG. NLPMIN3 - Nonlinear Programming user supplied grad. DN0ONF. NLSTART - Generate starting values for NL routines. QPMIN - Quadratic Programming. SOLVEFREE - Set frequency of freeing temp variables. Random Number Generation and Testing - IMSL Names preserved I_RNGET - Gets seed used in IMSL Random Number generators. I_RNSET - Sets seed used in IMSL Random Number generators. I_RNOPG - Gets the type of generator currently in use. I_RNOPT - Selects the type of uniform (0,1) generator. I_DRNSES - Initializes the table used by shuffled generators. I_DRNGES - Get the table used in the shuffled generators. I_DRNUN - Uniform (0,1) Generator I_DRNNOR - Random Normal Distribution I_RNBIN - Random integers from binomial distribution I_DRNGDA - Random integers from discrete distribution alias approach I_DRNGDT - Random integers from discrete using table lookup I_RNGEO - Random integers from Geometric distribution I_RNHYP - Random integers from Hypergeometric distribution. I_RNNBN - Negative binomial distribution I_DRNBET - Random numbers from beta distribution I_DRNCHI - Random numbers from Chi-squared distribution I_DRNCHY - Random numbers from Cauchy distribution I_DRNEXP - Random numbers from standard exponential I_DRNEXT - Random numbers from mixture of two exponential distributions I_DRNGAM - Random numbers from standard gamma distribution I_DRNGCT - Random numbers from general continuous distribution I_DRNLNL - Random numbers from lognormal distribution I_DRNNOA - Random normal numbers using acceptance/rejection I_DRNNOR - Random normal numbers using CDF method I_DRNSTA - Random numbers from stable distribution I_DRNTRI - Random numbers from triangular distribution I_DRNVMS - Random numbers from Von Mises distribution I_DRNWIB - Random numbers from Weibull distribution I_RNMTN - Random numbers from multinomial distribution I_DRNMVN - Random numbers from multivariate normal I_DRNSPH - Random numbers on the unit circle I_RNPER - Random perturbation of integers I_RNSRI - Index of random sample without replacement Screen I/O and Plot Subroutines CLS - Clear screen. FPLOT - Plot a Function GRAPH - High Resolution graph. GRAPHP - Multi-Pass Graphics Programing Capability GRCHARSET - Set Character Set for Graphics GRREPLAY - Graph replay and reformat command. OUTDOUBLE - Display a Real*8 value at a x, y on screen. OUTINTEGER - Display an Integer*4 value at a x, y on screen. OUTSTRING - Display a string value at a x, y point on screen. PLOT - Line-Printer Graphics SCREENCLOSE - Turn off Display Manager SCREENOPEN - Turn on Display Manager SCREENOUTOFF - Turn screen output off. SCREENOUTON - Turn screen output on. Statistical Functions BETAPROB - Calculate a beta probability. BINDF - Evaluate Binomial Distribution Function BINPR - Evaluate Binomial Probability Function BJ_IDEN - ACF and PACF Listing and or Plotting BLUS - BLUS Residual Analysis CHISQPROB - Calculate chi-square probability. DERF - Error function of real*8/real*16 argument. DERFC - Inverse of error function. FPROB - Probability of F distribution. HYPDF - Evaluate Hypergeometric Distribution Function HYPPR - Evaluate Hypergeometric Probability Function INVBETA - Inverse beta distribution. INVCHISQ - Inverse Chi-square distribution. INVFDIS - Inverse F distribution. INVTDIS - Inverse t distribution. MEAN - Average of a 1d object.vector. MEDIAN - Median of a real*8 object. NCCHISQ - Non central chi-square probability. NORMDEN - Normal density. POIDF - Evaluate Poisson Distribution Function POIPR - Evaluate Poisson Probability Function PROBIT - Inverse normal distribution. PROBNORM - Probability of normal distribution. PROBNORM2 - Bivariate probability of Nornal distribution. Q1 - Q1 of a real*8 object Q3 - Q3 of a real*8 object TDEN - t distribution density. TPROB - t distribution probability. VARIANCE - Variance of an object. System Subroutines DODOS - Execute a command string if under dos/windows. DOUNIX - Execute a command string if under unix. SYSTEM - Issue a system command. Time Series Functions ACF - Calculate autocorrelation function of a 1d object. CCF - Calculate the cross correlation function on two objects. DIF - Difference a series. FRACDIF - Fractional Differencing. GENARMA - Generate an ARMA series given parameters. MAKEJUL - Make a Julian date from a time series Variable Precision Math Subroutines and Functions Subroutines VPASET - Set Variable Precision Math Options Functions VPA - Variable Precision Math calculation Listing of Alias function names Name Replaced by LOG DLOG LN DLOG LOG10 DLOG10 EXP DEXP MOD DMOD MAX DMAX MAX1 DMAX1 MIN DMIN MIN1 DMIN1 INVERSE INV SIN DSIN COS DCOS R4TOR8 DFLOAT R8TOR4 FLOAT SQRT DSQRT SINH DSINH GAMMA DGAMMA COSH DCOSH CONJ DCONJ ATAN DATAN ATAN2 DATAN2 ARSIN DARSIN ARCOS DARCOS ABS DABS R8TOR4 SNGL R4TOR8 DBLE Index of Toolkit Routines contained in matrix2.mac Index of programs, subroutines and Functions ACF_PLOT - Simple ACF Plot AUTOCOV - Autocovariance BPF - Baxter - King MA Filter BPFM - Baxter - King MA Filter with missing data CFREQ - Determine Cumulative Frequency Distribution COINT2 - Cointegration Tests of Two Series COINT2LM - Cointegration Tests of Two Series, OLS, L1, MM COINT2M - Moving Cointegration of Two Series COINT2ME - Moving Cointegration of Two Series - Extended Args. COINT2M2 - Moving Cointegration of Two Series OLS, L1, MM COINT3 - Cointegration Tests of Three Series COINT3ME - Moving Cointegration of Three Series DATA_ACF - Calculate ACF and PACF Plots DATA2ACF - Calculate ACF and PACF Plots added argument DATAVIEW - View a Series Under Menu Control DO_SPEC - Display Periodogram and Spectrum DO2SPEC - Display Periodogram and Spectrum added argument DUD - Derivative Free Nonlinear Estimation FDIFINFO - Fractional Differencing Information FILTER - High Pass / Low Pass Filter using Real FFT FILTERC - High Pass / Low Pass Filter using Complex FFT FORPLOT - Forecast Plot using GRAPHP GARCH2P - Two Pass GARCH Using ARMA Command GARCH2PF - Two pass GARCH using ARMA Command with forecasting GTEST - Tests output of a ARCH/GARCH Model GWRITE - Save Objects in GAUSS Format using one file GWRITE2 - Save Objects in GAUSS Format using two files HP_BP_1 - Baxter-King & Hodrick-Prescott Filtering HP_BP_2 - Baxter-King & Hodrick-Prescott Filtering Moving Window HP_2 - Hodrick & Prescott Filtering of a Moving Window IRF - Impulse Response Functions of VAR Model LMTEST - Engle (1982) test for ARCH for a range of lags MARSPLOT - Automatically plot MARS Curves and Surface Plots MCLEODLI - McLeod-Li (1983) Linearity test MARQ - Estimation of a Nonlinear Model using Derivatives MINIMAX - Estimate MINIMAX with MAXF2 MISSPLOT - Plot of a series with Missing Data MQSTAT - Multivariate Q Statistic MOVEAVE - Moving average of a vector MOVEBJ - Moving Arima Forecast using AUTOBJ MOVEH82 - Moving Hinich 82 test MOVEH96 - Moving Hinich 96 test MOVEOLS - Moving OLS Calculation MOVEVAR - Moving Variance NLVARCOV - Calculates NLLSQ Variance Covariance OLSPLOT - Plot of Fitted and Actual Data & Res PAD - Pad a 1D Real*8 Series on both ends PVALUE_1 - Present value of $1 recieved at end of n years PVALUE_2 - Present Value of an Annuity of $1 PVALUE_3 - Present value of $1 recieved throughout year RTEST - Test Residuals of Model RTEST2 - Test Residuals of Model - No RES and Y Plots QUANTREG - Quantile Regression Program RESET77 - Thursby - Schmidt Regression Specification Test SUBSET - Subset 1d, 2d array, vector or matrix under a mask. Running the Matrix Command: The MATRIX command can be run in "batch mode" or interactively. When run interactively, statements can be reloaded and scripts can be edited and submitted. The file _imatrix.mac is the default name for the script files. The first member _matrix is automatically run if the script is submitted. The B34S save file format is the same as used by Speakeasy and the Speakeasy importall statement can be used to further process this file. The MATRIX command can also read and write SCA and RATS portable files and provides a bridge between these systems. The Display Manager provides a number of options on how to use the system. 1. The MATRIX button on the Display Manager window gets B34S into interactive MATRIX command mode. If the user wants to customize how the matrix command will come up, the IMATRIX macro in the MATRIX.MAC file can be changed. The user must not make an error in this file since the whole interative matrix command will be unable to execute if there is an error.. In MANUAL mode one command at a time can be entered or scripts can be edited and run. Interative use is intended for quick commands. If the IMATRIX member of the MATRIX.MAC file is modified, the MATRIX button could be made to execute a matrix program, or another B34S command BEFORE going into MANUAL mode in the matrix command. The default IMATRIX member of matrix.mac is: /$ This shell can be modified to load data set if desired b34sexec matrix; call manual; b34srun; In the user directory it is possible is the _imatrix.mac file that has a required member _matrix which can be edited and submitted while under in the interactive matrix command mode. A sample file is: ==_matrix program _matrix; /$ Lines above this are required x=array(100:); y=rn(x); z=rn(y); x=rn(x); call graph(y,z); call names(all); call tabulate(x,y,x); call olsq(y x:print); call load(user '_imatrix.mac'); call user; /$ Lines next are required call manual; return; end; == ==user program user; /$ this is a user program; call print('I am in the user program user!!':); return; end; == A simple example is: /; /; Sample user _imatrix.mac file /; ==_MATRIX program _matrix; /; Commands after here x=rn(matrix(5,5:)); ix=inv(x); call print(x,ix,x*ix); call print('How are things going'); return; end; == Remember the _imatrix.mac file is in the user directory. It submits matrix program while inside the matrix command. The imatrix macro in matrix.mac allows data to be loaded etc BEFORE calling the matrix command!! 2. The MATRIX command has been designed to run in "batch" mode from a command file. This file is usually submitted from a FILE or from the TASKS command of the Display Manager. 3. The LINEEDIT option under the MENU button allows multiple lines to be entered. This way quick MATRIX jobs can be submitted. By use of the MENU "RELOAD" facility, the commands just submitted can be recalled and changed but are not saved in memory. 4. The MENU command provides two options. The MATRIX command gives the user a screen where any MATRIX commands can be entered. In this mode a "batch file" is interactively submitted. This mode is designed for quick jobs that call user PROGRAMS, SUBROUTINES and FUNCTIONS. By use of the RELOAD option, these commands can be modified and resubmitted. The MENU command also allows access to the MATRIX command which is the same as what is available with the IMATRIX Display Manager command (see # 1 above). In the middle of a MATRIX command section the command: call save; will save the workspace. Later in the session OR at a later session, the command: call restore; will bring back what has been saved. The default name is matrix.psv A specific file can be given using the forms call save(:file 'mystuff.psv'); which can be brought back with call restore(:file 'mystuff.psv'); These save files can be read with Speakeasy(r) using the Speakeasy IMPORTALL command. If b34s matrix command PROGRAMS, SUBROUTINES or FUNCTIONS are present in the B34S MATRIX command workspace, the command call save(:speakeasy); will save only data objects. Because Speakeasy will not support variables of the form %name, such variables are name _name. On a restore these are converted back. The keyword checkpoint can be used in the place of save to be compatible with Speakeasy. If checkpoint is used, subroutines, functions and programs are automatically not saved. The MATRIX command interactive mode can be terminated by call run; and started later in the same job with call manual; This way a complex program can be debuged by looking an intermediate values. Since all output is written to the b34s output file in the usual manner, View log and View output are always available to inspect results. Help The usual b34s help facility provides a discussion of all B34S MATRIX commands. In addition the file matrix.mac provides an example of virtually all MATRIX commands. These examples can be viewed with the help facility, executed or saved in the B34S buffer where they can be either run or further modified. The file matrix2.mac contains programs, subroutines or fuctions that can be run provided that they are loaded with the command call load(somename); Form of MATRIX command. b34sexec matrix options parameters; b34seend$ MATRIX Command Options SAVEASVECTOR - Saves B34S variables as vectors. Default is to save as an array. For the differences between vector and array math, see below. For a minor example, if the b34s variable GASOUT is saved as an array. The statement test=gasout**2.; will square each element, If GASOUT is saved as a vector, the command test=gasout**2.; will produce one observation which is the sum of squared values of GASOUT. SHOWUSE - Shows use of arrays. Indicates current limits of the program. HEADER - Turns on Header. Default is NOHEADER. The commands call noheader; call header; can turn the header off and on inside the matrix step. CBUFFER=n1 - Sets size of Command buffer. If a bigger size is needed, a message will be given. For further info see below. SIGD=n2 - Sets number of significant digits in variable printing. Default g16.8 for real*8. This can also be set by call sigd(i); in open matrix language code. If i le 8 g16.i is used. If i > 8 then g16+(i-8).i format is used. For example i = 10 uses g18.10. SOLVEFREE=n3 - Sets frequency of cleaning temp variables with SOLVE command. Default = 50. If error message on exceeding temp variables is seen, set number smaller. A larger n2 value can reduce CPU time. The subroutine SOLVEFREE can be used to test your code. DSEED=r - Sets the seed for GGNML and GGUBS and other routines. This seed works the same but can be set to a different value than the seed set under the OPTIONS command. The exact random number routine used is set on the OPTIONS command with the command RECVER and RNVER. Consult help files for this command for further detail. Default = 123457.0d+00 or the seed set with SETSEED on the OPTIONS command. The possibly of a different seed for the MATRIX command isolates how this command runs in relation to a default seed set with SETSEED in a user AUTOEXEC.B34 file. Since not all random number generators are equal, users should use caution. Empirically it was found that if 3000 random normal deviates are generated using the default (IMSL) seed, the BDS statistic flags the series as "nonlinear." This finding is strange to say the least. Using I_RNOPT command the exact uniform generator for IMSL random number routines can be selected. The MATRIX command RN( ) and REC( ) have optional switches that allow the default random number routines to be replaced by IMSL version 10 routines. DISPLAY=key - Sets accuracy and number of Cols to print matrix output. Key can be set: COL80FIXED f18.5 COL80MEDIUM g15.6 COL80HIGH g25.16 COL129FIXED f18.5 COL129MEDIUM g15.6 COL129HIGH g25.16 The default is COL129MEDIUM. This can be set globally with the options commands linesize and sigd LINESIZE SIGD DISPLAY 80 < 8 COL80FIXED 80 8 COL80MEDIUM 80 > 8 COL80HIGH 132 < 8 COL129FIXED 132 8 COL129MEDIUM 132 > 8 COL129HIGH ADJUST_PARSE - Can be used if functions are being called outside subroutines only. For more detail see the FUNCTION help file. This option is not the the general user. Size of problems: The maximum size of any problem is limited by the maximum number of sentences in the parse table (MAXNSENT) and the maximum number of tokens (MAXNTOKEN). These can be increased in the autoexec.b34 file on the PC, or in an OPTIONS command such as b34sexec options maxnsent(600) maxntoken(5000); b34srun; Other current constraints are: Max number of objects in memory 10000. Max number of arguments to a function or subroutine 2500. These constraints may change in the future. For other limitations see the SHOWUSE command of the MATRIX sentence. Overview of the Matrix Language Cont. The goal of the MATRIX command is to provide a high-level programing language that will allow the user to easily customize a B34S application. The MATRIX command consists of analytic statements that that can be reduced to assingment statements and keywords. In contrast to address oriented programing languages such as FORTRAN, the B34S matrix language is object oriented. Assuming X was a vector of 20 numbers, the command ameanx=mean(x); places the mean of the x vector in the variable ameanx. Keywords such as PRINT are used in the form call print('This will print x',x); which will print the statement and the object x. The PRINT command will allow up to 400 arguments or objects to be passed in one call. Keywords are built into the language or can be user written SUBROUTINES, FUNCTIONS or PROGRAMS. User SUBROUTINES and FUNCTIONS pass objects as arguments and calculate in a protected workspace. This means that a variable X can be defined at both the base level and in the SUBROUTINE to mean different things. The difference between PROGRAMS and SUBROUTINES is that SUBROUTINES pass arguments and work in a protected space while PROGRAMS do not pass arguments. PROGRAMS use the current level as their name domain. Hence if a PROGRAM is called from the base level, it can access objects known at that level. If the same programn is called from a SUBROUTINE, it will access data known at that level only. User SUBROUTINES and FUNCTIONS also have access to GLOBAL variables. Recursive Notes: User SUBROUTINES and FUNCTIONS can be called recursively up to the MAXSTAK limit currently 5000, although CBUFFER will have to be increased to handle saving the open arguments. Recursive calls are very very slow but are useful in some situations. Job RECURSIVE in matrix.mac shows that for the same task, a DO loop is 100 times faster than a recursive function or subroutine call. A DO loop itself is slow relative to object calculations. FORMULA & SOLVE statements are from 4 - 10 times faster than DO loops. The speed gain depends on: 1. the length of the vectors, 2: the complexity of the problem and 3. the SOLVEFREE setting which sets how frequently temp variables are freed in the workspace. If recursive calls are a major problem, a brach to an externallu compiled Fortran program might be a good solution. For further detail, see section 14 of this help file. User SUBROUTINES, PROGRAMS and FUNCTIONS must have names of 8 or less characters. A simple example of a user subroutine is DESC subroutine desc(name,mean1,var); mean1=mean(name); var=variance(name); return; end; The following code fragment uses DESC x=array(5:1 2 3 4 5); call desc(x,m,v); call print(Variable mean var',x,m,v); A PROGRAM for the same result would be program desc2; m=mean(x); v=variance(x); return; end; The following code fragment uses DESC2. x=array(5:1 2 3 4 5); call desc2; call print(Variable mean var',x,m,v); Note that in this case the actual names are needed to be used inside DESC2 since the name domain is at the global level. A simple user function DESF will return the mean although in this case it makes little sense since we have the built-in function mean. function desf(x); m=mean(x); return(m); end; The following code fragment uses desf x=array(5:1 2 3 4 5); amean=desf(x); MATRIX command built-in FUNCTIONS and SUBROUTINES in addition to allowing name arguments, also allow passing character strings subject to a limitation shown below. Given subroutine printit(c); call print(c); return; end; If the string is LE 8 characters the form call printit('less8'); can be used. If the string is greater that 8, then the correct form is call character(c,'This is greater than 8'); call printit(c); must be used. This is illustrated in the job char_4 in the matrix.mac file. What is happening is that the form call printit('aa'); places the string 'aa' in a character*8 variable. The form call character(cc,'aa'); call printit(cc); places the string 'aa' in a two element character*1 array. Thus the statement name='J'; creates a character*8 variable. If this is used in an assingment statement to a character*1 variable, the character*1 variable will be redefined which may not be what is wanted. This is illustrated in the job char_5 which is listed next: b34sexec matrix; call echooff; x=rn(matrix(30,30:)); call character(ii,'Element (1, '); call character(ii2,' '); jj=integers(12,17); do i=1,30; call inttostr(i,ii2,'(i6)'); ii(jj)=ii2(jj-11); call character(rp,')'); /$ ********************************************************** /$ Warning the statement /$ ii(18)=')'; /$ does not work since it will be redefined to be character*8 /$ and will be outside the 132 range and not printed /$ ********************************************************** ii(18)=rp; call print(ii,x(1,i) :line); next i; call names(all); b34srun; Summary. Since the MATRIX command supports character*8 and character*1 data added care must be used to pass just what is desired to a user routine. Built-in subroutines and functions can detect just what has been passed and take the appropriate path. The command character or char1 can be used to build a multi-line 2D character*1 array. call character(text,'This is line one' 'This is line two which is longer' 'This is line three' 'This is 4'); produces a 2d character*1 array of size 4,132 where 132 is the max length of a line. Strings to the character command can be up to 132 in length. Advanced concepts. Note: If a MATRIX command variable name is passed to a user SUBROUTINE or FUNCTION, it is possible to return a changed value. If a structured object is passed, it is NOT possible to return a value. This follows the Speakeasy convention but is not the usual case in Fortran. The reason for this is that the structured object may repackage the data, which is not the case in Fortran. For example assume matrix A defined as n=10; a=rn(matrix(n,n:)); if we call desc (see routine listed above) as call desc(a(1,),m,v); mm(1)=m; vv(1)=v; call desc(a(2,),m,v); mm(2)=m; vv(2)=v; we have the desired result since DESC does not change the first argument. However the code call desc(a(1,),m(1),v(1)); call desc(a(2,),m(2),v(2)); will not work as intended because in the parsing process, m(1), m(2), v(1), and v(2) are replaced by temporary variables. After returning the temp values not the real values are saved. One would think that this could be taken care of by the compiler BUT this is not possible. The reason relates to how stuctured objects are handled. Consider the following code. x=rn(matrix(10,10:)); row2=x(2,); In the matrix the elements of row 2 differ by 10 memory positions since Fortran saves by column. However the command row2=x(2,); "repackages" these elements so that they are next to each other. Hence putting them back into the matrix X is not just a copy operation. For this reason users have to be careful when passing structured objects. B34S MATRIX USER SUBROUTINES and FUNCTIONS pass structured objects by value NOT by address! More detail on this is provided below. In other languages such as MATLAB arguments cannot be changed at all. In MATLAB if the user has a function myfun that takes three arguments and returns two, the MATLAB calling form is [back1,back2]=myfun(in1,in2,in3); where in matlab the function is defined by function[bb,cc]=myfun(aa1,aa2,aa3) The B34S MATRIX command follows the Speakeasy conventions and allows a bit more flexibility. However users have to be careful. The MATRIX language recognizes the following types of objects which are given kind and klass numbers. kind character*1 -1 real*4 4 real*8 8 real*16 -16 complex*16 16 complex*32 32 program 1 subroutine 2 function 3 formula 33 character*8 -8 integer*4 -4 integer*8 -88 vpa real 88 vpa real packed 888 vpa integer -44 vpa integer packed -444 vpa complex 160 vpa complex packed 1600 not defined -99 Klass scalar 0 vector 1 matrix 2 1 dim array 5 2 dim array 6 The matrix language recognizes structured objects. Define a as matrix 1. 2. 3. 4. 5. 6. 7. 8. 9. with the command a=matrix(3,3:1 2 3 4 5 6 7 8 9); The matrix, array and vector commands are the only ones that automatically convert integer input to real*8. If an integer*4 matrix is desired use a=idint(matrix(3,3:1 2 3 4 5 6 7 8 9)); The command call print(a(2,)); prints row 2 while call print(a(,3)) prints column 3. After defining an integer vector i of two elements 1 and 3 i=integers(1,3,2); the commands call print('This is row 1 and 3',a(i,)); call print('This is col 1 and 3',a(,i)); will pull off rows 1 & 3 and columns 1 & 3 respectively. Remember structured objects pass as values NOT addresses. Structured objects can be used as arguments but cannot be used for output. On the left of an = sign assingments are possible. This is contrary to Fortran which passes by address. B34S "repackages" the structured object data and passes the values leaving the original array intact. Internally the matrix is saved by columns following Fortran. Hence both x(2,) and x(,2) pass contiguous values although the underlying storage is NOT contiguous for x(2,). What is happening is that the B34S parser "repackages" the second row, which was not saved as contiguous values, into a vector. For this reason it is not possible to return back values from a structured variable address such as x(,i), x(i) or x(3,) etc. In a sense this "limitation" of returning values in a structured object frees the user from having to remember the underlying storage of the matrix object. Here the vector is an object in its own right. That is how we think of it in mathmatics. The goal of the MATRIX design is to think in mathematics, not in computer storage conventions. Speed design considerations: Code such as i=integers(2,20); j=i-1; a(j) = b(i); involves repetitive expansion of A which has not been defined to the correct size and therefore has to be expanded on the fly. While such code will work, a better approach is b34sexec matrix; b=rn(array(30:)); i=integers(2,20); j=i-1; a=array(dmax(j):); a(j) = b(i); call tabulate(a,b); b34seend$ which will run substantially faster. The following code illustrates some of these concepts. x=rn(matrix(3,3:)); a(1,1)= mean(x); v=vector(3:1 2 3); a(2,)=v; i=idint(array(2:1,3)); * place terms 1 3 of v in newv; newv=v(i); x=matrix(3,3:integers(1,9)); * places rows 1 and 3 of x or 1 2 3 7 8 9 in newx; newx=x(i,); The mean of x is now in a(1,1), 1 2 3 is in row 1 of newx and 7 8 9 are in row 2. 1.0 and 3.0 are in variable newv. Data can be placed in structured objects a number of ways. To put 6.0 in col 3 of x at all locations: call setcol(x,3,6.0); To put 7.0 in row 6 of x at all locations: call setrow(x,6,7.0); another way would be x(6,)=(array(nocols(x):)+7.0); assuming x is a 2d array. Or x(6,)=7.; Object expansion is possible. The command a(6,)=v; would place v in a new row 6 of x. Due to possible ambiguity, if v is a vector and i and j are vectors, the statement x(i,j)=v; is not allowed. The command newx=matrix(n,n:)+8.0; will build a n by n matrix with 8.0 on the diagonal while newx=array(n,n:)+8.0; will put 8.0 in all elements. These statements illustrate the difference between "matrix math" and "array math." The operators .EQ. .LT. .GT. .NE. .GE. .LE. .OR. .AND. can be used in expressions as well as IF and WHILE statements. Given x=2.0; y=3.0; test1=x.eq.y; test2=x.ne.y; returns test1=0.0 and test2=1.0 If x=array(:1 2 3); y=array(:1 2 4); the commands test1=x.eq.y; test2=x.ne.y; imply test1=array(:1.0 1.0 0.0) since x(3) ne y(3) test2=array(:0.0 0.0 1.0) since x(1)=y(1) and x(2)=y(2) Note: That in an IF statement the arguments must be scalars while analytic statements of the form x= will create a scalar or vector depending on the inputs. If one input is a vector, the other input must be a scalar or a vector of the same size. Warning: The logical operators .and. and .or. can be used to compare 0.0 and 1.0 values only. Names of Objects and commands. B34S MATRIX commands that do not return an argument across an equals are executed by the CALL sentence. The CALL sentence first looks in named storage for a routine with this name. If this is not found, then the built in routines are used. While it is possible to have a user routine with the same name as a built in routine, this is not a good idea. For example assume you have loaded the series GASOUT. The command acf=acf(gasout,24); will place 24 autocorrelations of the GASOUT series in the structured variable ACF. The result is that the ACF command is no longer available. If the command acf2=acf(gasin,24); is given later in the job, B34S will object that ACF is not a 2D object. The solution is to use FREE to "turn on" the built-in command ACF. call free(acf); acf2=acf(gasin,24); or better still resist using an internal command name as a variable name. Math using the MATRIX Command. The MATRIX command allows math between 1 and 2 dimensional arrays and between vectors (1 dimensional) and matrices (2 dimensional) for real*8 and complex*16 objects. Array math can be done between real*8, real*4, real*16, integer*4 and complex*16 objects. Objects of different types cannot be mixed. If this is attempted, B34S will give a "mixed mode" error message. The reason for this is given below. Warning. Only in the most extreme cases should variables be saved as real*4. The real*4 data type is for some graphic applications and for variable storage in cases where memory is low. While array and matrix math is possible for real*4, other functions such as sin will not work to prevent users from making calcualtios that have accyracy loss. Matrix and vector math can be done only with real*4 real*8, real*16, complex*16, complex*32 variables and vpa real and complex objects. Assume r8 is real*8. newr=r8*2.0; is allowed but newr=r8*2; is not allowed since 2 is an integer. The error message will refer to "mixed mode" operations. The logic behind not allowing r8*2 is that it is not clear whether an integer*4 or real*8 result is desired. This convention is not followed in Speakeasy which tries to make everything real*8. The reason for this "restriction" is illustrated by the following example. Assume the commands x=10.; y=x*2; where given. The parser would not know if Y should be integer*4 20 or real*8 20. Note: The functions DBLE, SNGL, DINT, IDINT, INT, FLOAT, DFLOAT, REAL, IMAG, COMPLEX, r8tor16, r16tor8, c16toc32, c32toc16 i4toi8, i8toi4 and vpa can be used to convert the storage of a variable. Speed considerations: To increase run-time speed the CALL command is used in place of just naming the SUBROUTINE as is done with some other programming systems such as Speakeasy and MATLAB. B34S follows the Fortran convention. Arguments passed to subroutines and functions must be inside ( ). Built in Matrix Commands called with CALL must not be mixed. For example call names; is allowed. but call names names(all); is not allowed since the command NAMES(ALL) is not the first key word in the call sentence. Matrix Command Files: The file c:\b34slm\matrix.mac contains a number of complete matrix command examples. These can be executed, loaded into the program buffer under TASKS, or viewed using the help facility. There are working examples for virtually every matrix command. Users should customize these examples. The file c:\b34slm\matrix2.mac contains matrix command SUBROUTINES, PROGRAMS and FUNCTIONS. These can be loaded as run time with the command call load(somename); or loaded into the program buffer, or viewed with the help facility. Each file cannot contain anything except, SUBROUTINES, PROGRAMS or FUNCTIONS. The advantage of loading a routine with a CALL LOAD statement, rather than having the routine in the command file, is that parse space is saved and execution speed is also increased. If the user adds routines to the library, it is important that the code be "clean" since the parser will not check the code in the same manner as would be the case if the routine were loaded. Unless the routine produces output, most routines are called inside CALL ECHOOFF; CALL ECHOON; statements so that the statements executed will not echo to the output file. The files staging.mac and staging2.mac are like matrix.mac and matrix2 except that they are for prospective commands. These commands are documented inside the routine and not in the b34shelp.dat file. Basic rules of the B34S Matrix Command: 1. All matrix statements MUST end in $ or ;. x=dsin(q); 2. Mixed mode math is not allowed. For example assuming x is real*8 x=x*2; is not allowed because x is real*8 and 2 is an integer. The reason mixed mode is not allowed is that the processor would not know what to do with the result. The correct form is x=x*2.; or x=idint(x)*2; if you want an integer result and x was real*8 before the command. 3. Calculated Structured objects can only be used on the right of an expression or in a subroutine call as input. For example mm=mean(x(,3)); calculates the mean of col 3 while nn=mean(x(3,)); calculates the mean of row 3. i=integers(2,30); y(i)=x(i-1); copies x(1),...,x(29) to y(2),...,y(30). The code below is not allowed since it involves a calculated left hand side i=integers(1,29); y(i-1)=x(i); 4. Structured objects can be used on the left of an assignment statement to load data. The commands x=3.0; x(2)=4.0; add another element to x. x=rn(matrix(4,4:)); x(,2)=0.; places 0.0 in col 2 while x(3,)=99.; places 99.0 in row 3. The following code shows advanced structured index processing. This code is available in matrix.mac as overview_2 /$ Illustrates Structural Index Processing b34sexec matrix; x =rn(matrix(6,6:)); y =matrix(6,6:); yy =matrix(6,6:); z =matrix(6,6:); zz =matrix(6,6:); i=integers(4,6); j=integers(1,3); xhold=x; hold=x(,i); call print('cols 4-6 x go to hold',x,hold); y(i, )=xhold(j,); call print('Rows 1-3 xhold in rows 4-6 y ',xhold,y); y=y*0.0; j2 =xhold(j,); y(i, )=j2 ; call print('Rows 1-3 xhold in rows 4-6 y ',xhold,y); z(,i)=xhold(,j); call print('cols 1-3 xhold in cols 4-6 z ',xhold,z); j55 =xhold(,j); z=z*0.0; z(,i)=j55; call print('cols 1-3 xhold in cols 4-6 z ',xhold,z); yy=yy*0.0; yy(i,)=xhold; call print('rows 1-3 xhold in rows 4-6 yy',xhold,yy); zz=zz*0.0; do ii=1,3; jj=ii+3; zz(,jj)=xhold(ii,); enddo; /; i=integers(4,6); /; j=integers(1,3); call print('Note that zz(,i) = xhold(j,) will not work':); call print('Testing zzalt(,i)= transpose(xhold(j,))':); /; Use of Transpose speeds things up over do loop zzalt=zz*0.0; zzalt(,i)= transpose(xhold(j,)) ; call print('rows 1-3 xhold in cols 4-6 zz',xhold,zz,zzalt); zz=zz*0.0; zzalt=zz; do ii=1,3; jj=ii+3; zz(jj,)=xhold(,ii); enddo; call print('Note that zz(i,) =xhold(,j) will not work':); call print('Testing zzalt(i,)= transpose(xhold(,j))':); zzalt(i,)=transpose(xhold(,j)); call print('cols 1-3 xhold in rows 4-6 zz',xhold,zz,zzalt); oldx=rn(matrix(20,6:)); newx= matrix(20,5:); i=integers(4); newx(,i)=oldx(,i); call print('Col 1-4 in oldx goes to newx',oldx,newx); oldx=rn(matrix(20,6:)); newx= matrix(20,5:); i=integers(4); newx(1,i)=oldx(1,i); call print('This puts the first element in col ',oldx,newx); newx=newx*0.0; newx(i,1)=oldx(i,1); call print('This puts the first element in row ',oldx,newx); newx=newx*0.0; newx( ,i)=oldx( ,i); call print('Whole col copied here',oldx,newx); oldx=rn(matrix(10,5:)); newx= matrix(20,5:); i=integers(4); newx(i,1)=oldx(i,1); call print('This puts the first element in row ',oldx,newx); newx=newx*0.0; newx(i,)=oldx(i,); call print('Whole row copied',oldx,newx); * We subset a matrix here ; a=rn(matrix(10,5:)); call print('Pull off rows 1-3, cols 2-4', a,a(integers(1,3),integers(2,4))); b34srun; 5. Functions or math expressions are not allowed on the left hand side of an equation. Assume the user wants to load another row. The command x(norows(x)+1,)=v; in the sequence x=matrix(3,3:1 2 3 4 5 6 7 8 9); v=vector(3:22 33 44); x(norows(x)+1,)=v; will not work. The correct way to proceed is: x=matrix(3,3:1 2 3 4 5 6 7 8 9); v=vector(3:22 33 44); n=norows(x)+1; x(n,)=v; to produce 1. 2. 3. 4. 5. 6. 7. 8. 9. 22. 33. 44. The command x(i+1)=value; will not work since there is a calculation implicit on the left. The correct code is: j=i+1; x(j)=value; 6. Matrix and array math is supported. If x is a 3 by 3 matrix, the command x=afam(x); will create the 3 by 3 array x. If x is a 3 by 1 array. The command mx=vfam(x); will create a matrix mx containing columns of x. If x=rn(array(3,3:)); vv=vfam(x); vv is a 3 by 3 matrix. To convert x to a vector column by column use vvnew=vector(:x); 7. Keywords should not be used as variable names. If they are the command with this name is "turned off." This can cause unpredictable results with user PROGRAMS, SUBROUTINES and FUNCTIONS. Keywords cannot conflict with user program, subroutine or function names since the users code is not loaded unless a statement of the form call load(name); is given. 8. MATRIX command SUBROUTINES and FUNCTIONS allow passing arguments. For example: call myprog(x,y); y=myfunc(x,y); Character values can be passed and optionally changed. b34sexec matrix; subroutine test(a); call print('In routine test A= ',a); call character(a,'This is a very long string'); return; end; /$ pass in character*8 call test('junk'); call character(jj,'some junk going in'); call print(jj); /$ pass in a character*1 array call test(jj); call print(jj); b34srun; Special characters such as : and | are not allowed in USER SUBROUTINES or FUNCTION calls because of the difficulty of parsing these characters in the user routine. This restriction may change in future versions of the MATRIX command if there is demand. 9. Coding assumptions. Statements such as: y = x-z; are allowed. Statements such as y = -x+z; will not work as intended. The error message will be "Cannot classify sentence Y ...". The command should be given as y = -1.*x + z; or better still y = (-1.*x) + z; A statement y = x**2; where x is real*8 will get a mixed mode message and should be given as y = x**2.; Complex statements such as yhat = b1*dexp(-b2*x)+ b3*dexp(-(x-b4)**2./b5**2.) + b6*dexp(-(x-b7)**2./b8**2.); will not work and should have ( ) around the power expressions and -1.* . yhat = b1*dexp(-1.0*b2*x)+ b3*dexp(-1.0*((x-b4)**2.)/(b5**2.)) + b6*dexp(-1.0*((x-b7)**2.)/(b8**2.)); Examples of MATRIX language statements: The statement y=dsin(x); is an analytic statement that creates the stuctured object y by taking the sin of structured object x. The statements below defines a user subroutine moveaver to calculate a moving average. subroutine moveaver(x,nterms,moveaver); 10. The following code illustrates automatic expansion. x(1)=10.; x(2)=20.; The array x contains elements 10. and 20. Warning: The commands x(1)=10.; x(2)=20; produces an array of 0 20 since the statement x(2)=20 redefines the x array to be integer! This is an easy mistake to make! Computers do what we tell them to do! Statements such as x(0) = 20.; x(-1)= 20.; x(1) = 20.; all set element 1 of x to 20. The x(0) and x(-1) statements will generate a message warning the user. 11. Memory Management Issues. Warning: Automatic expansion inside a subroutine, program, DO loop, dowhile loop, or function can cause the program to "waste" memory since newer copies of the expanded variable will not fit into the old location. The matrix command will have to allocate a new space which will leave a "hole" in memory. B34S provides the capability of the user programming memory management as needed. The command call compress; can be used to compress the workspace in this situation. A better solution is to allocate the variable to the max shape outside the loop. In addition to space requirements, prior allocation will substantially speed up execution. If memory problems are encountered, the command call names(all); can be used to see how the variables are saved in memory and whether as the loop proceeds more space is used. In some situations the max temp number message may be reached if temps cannot be automatically freed. For example compare the following code; n=10; x=array(n:); call names(all); do i=1,n; x(i)=dfloat(i); call names(all); enddo; with n=10; call names(all); do i=1,n; x(i)=dfloat(i); call names(all); enddo; The first job will run faster and not use up memory. This job can be found in the file matrix.mac under MEMORY and should be run by users wanting to write efficient subroutines. The command n=10; call compress(n); allows the user to compress every n times call compress is called. For example: do i=1,nn; * statements here; call compress(100); enddo; Compresses every 100 times the loop calls compress. An alternative is to use the solvefree command as do i=1,2000; call solvefree(:alttemp); * many commands here ; call solvefree(:cleantemp); enddo; Due to enhancements to the memory management in B34S that were part of the 8.10G version, this option only has use if memory management is needed as part of a user model that is called by one of the nonlinear commands. For added detail on this use, see below. In the above example the first call with :alttemp sets %%____ style temp variables in place of the default ##____ style. The command :cleantemp resets the temp style to ##____ and cleans all %%____ temps, leaving the ##_____ style temps in place. If this capability is used carefully, substantial speed gains can be made. In addition the max number of temps will not be reached. Use of this feature slows down processing and is usually not needed. The command call solvefree(:cleantemp2); cleans user temps at or above the current level. This can be useful within a nested call to clean work space. The dowhile loop usually is cycled many times and may need active memory management. An example is: b34sexec matrix; sum=0.0; add=1.; ccount=1.; count=1.; tol=.1e-6; /$ outer dowhile does things 2 times call outstring(2,2,'We sum until we can add nothing!!'); call outstring(2,4,'Tol set as '); call outdouble(20,4,tol); call echooff; call solvefree(:alttemp); dowhile(ccount.ge.1..and.ccount.le.3.); sum=0.0; add=1.; count=1.; dowhile(add.gt.tol); oldsum=sum; sum=oldsum+((1./count)**3.); count=count+1.; call outdouble(2,6,add); add=sum-oldsum; /$ This section cleans temps if(dmod(count,10.).eq.0.)then; call solvefree(:cleantemp); call solvefree(:alttemp); endif; enddowhile; ccount=ccount+1.; call print('Sum was ',sum:); call print('Count was ',count); enddowhile; b34srun; Note: The command call compress; will be ignored if it is used in a program, function or subroutine that is called as part of a nonlinear estimation command such as NLLLSQ, CMAX2 etc. The reason for this restriction is to avoid the possibility of data movement that is not known to the calling command. If memory management is needed in this case, use the solvefree command. 12. Missing data Handling missing data often causes problems. Assume the following code: x=rn(array(10:)); lagx=lag(x,1); y=x-(10.*lagx); goody=goodrow(y); call tabulate(x,lagx,y,goody); Y will contain missing data in row 1. The variable goody will contain 9 non missing value observations. 13. Recursive solutions In many cases the solution to a problem requires recursive evaluation of an expression. While the use of recursive function calls is possible, it is not desirable since there is great overhead in calling the function or subroutine over and over again. The do loop, while still slow, is approximately 100 times faster than the recursive function call. The test problem RECURSIVE in c:\b34slm\matrix.mac documents how slow the recursive function call and do loop are for large problems. Another reason that a recursive function call is not recommended is that the stack must be saved. One way to handle a recursive call is to use the SOLVE statement to define the expression that has to be evaluated one observation at a time. If the expression contains multiple expressions that are the same, a FORMULA can be defined and used in the SOLVE statement. The FORMULA and SOLVE statements evaluate an expression over a range, ONE OBSERVATION AT A TIME. This in contrast to the usual analytic expression which is evaluated completely on the right BEFORE the copy is made. Unlike an expression, a formula or SOLVE statement can refer to itself on the right. The BLOCK keyword determines the order in which the formulas are evaluated. If the expression in the SOLVE statement does not have duplicate code, it is faster not to define a FORMULA. Examples of both approaches are given next. The code: test=array(10:); test(1)=.1; b=.9; solve(test=b*test(t-1)+rn(1.) :range 2 norows(test) :block test); call print(test); works but test = b*lag(test,1)+rn(1.); will not get the "correct" answer since the right hand side is built before the copy is done. More complex expressions. The FORMULA statement requires use of the subscript t unless the variable is a scalar. The use of the FORMULA and SOLVE statements are illustrated below: Example 1. b34sexec options ginclude('gas.b34'); b34srun; b34sexec matrix ; double=array(norows(gasout):); formula double = dlog(gasout(t)*2.); call names; call print(double); test2=array(norows(gasout):); solve(test2=test2(t-1)+double(t)+double(t-1) :range 2, norows(gasout) :block double); call print(mean(test2)); b34srun; Example 2. The following two statements are the same do i=1,n; x(i)=y(i)/2.; enddo; solve(x=y(t)/2. :range 1 n); Note: The SOLVE and FORMULA statements cannot use user functions. A DO loop can use user functions. More detail on the SOLVE and FORMULA statements are given below. An alternative to the FORMULA/SOLVE approach is to call a Fortran or other execuitable directly for the recursive calculation from inside a B34S Matrix command program. For an example of this approach see section 22 below and the examples FORTRAN and FORTRAN2 in the file matrix.mac. This method is very fast and is the best way to go. The down side is building the fortran and having to have a fortran compiler. Since c routines can also be used and there are GPL c compilers, this is not a serious problem, especially on linux where f77 cones with the operating system. 14. User defined data structures The B34S MATRIX command allows users to build custom data types. The below example shows the structure PEOPLE consisting of a name field (PNAME), a SSN field (SSN), an age field (AGE), a race field (RACE) and an income field (INCOME). The built in function SEXTRACT( ) is used to take out a field and the built in SUBROUTINE ISEXTRACT is used to place data in a structure. Both SEXTRACT and ISEXTRACT allow a third argument that operates on an element. The name SEXTRACT is "structure extract" while ISEXTRACT is "inverse structure extract." Use of these commands is illustrated below: b34sexec matrix; people=namelist(pname,ssn,age,race,income); pname =namelist(sue,joan,bob); ssn =array(:20,100,122); age =idint(array(:35,45,58)); race =namelist(hisp,white,black); income=array(:40000,35000,50000); call tabulate(pname,ssn,age,race,income); call print('This prints the age vector',sextract(people(3))); call print('Second person',sextract(people(1),2), sextract(people(3),2)); * make everyone a year older ; nage=age+1; call isextract(people(3),nage); call print(age); * make first person 77 years old; call isextract(people(3),77,1); call print(age); b34srun; Data structures are very powerful and, in the hands of an expert programmer, can be made to bring order to complex problems. 15. Advanced programming Concepts and Techniques for Large Problems Programs such as Speakeasy which are meant to be used interactively have automatic workspace compression. As a result a Speakeasy LINKULE programmer has to check for movement of defined objects anytime an object is created or freed. So as to increase speed, B34S was designed not to move the variables inside the allocator unless told to do so. If a DO loop terminates and the user is not in a SUBROUTINE, temp variables are freed. If a new temp variable is needed, B34S will try to place this variable in a free slot. If a variable is growing, this may not be possible. Hence it is good programming practice to create arrays and not rely on automatic variable expansion. In a subroutine call, a variable passed in is first copied to another location and set to the current level + 1. Thus there are name storage implications of a subroutine call. In a large program the command call compress; will manually clean out all temp variables and compress the workspace. While this command takes time, in a large job it may be required to save time and space. For example: temp variables are named ##1 ...... ##999999. If the peak number of temp variables gets > 999999, then B34S has to reuse old names and as a result slows down checking to see if a name is currently being used. A call to COMPRESS will reset the temp variable counter as well as free up space. Compacted space in worekingmemory speeds up execution. If call compress is called from a place it cannot run, say in a program subroutine or function that is being called by a b34s command such as cmax2, it will be ignored and no message will be given. The compress command will also be ignored if it is called under a running function, directly or indirectly. This means that a call to compress from a subroutine may or may not compress depending on whether the subroutine was called from a user function. The MATRIX command termination message gives space used, peak space used and peak and current temp # usage. Users can monitor their programs with these measures to optimize performance. In the opinion of the developer, the B34S MATRIX command DO loop is too slow. The problem is that the DO loop will start to run without knowing the ending point because it is supported at the lowest programming level. In contrast, Speakeasy requires that the user have a DO loop only in a SUBROUTINE, PROGRAM or FUNCTION where the loop end is known in theory. Ways to increase DO loop speed are high on the "to do" list. Faster CPU's may be the answer. The LF95 compiler appears to make faster DO loops than the older LF90 compiler. This suggests that the cache may be part of what is slowing things down. The test problem SOLVE6 illustrates some of these issues. Times and gains from the SOLVE statement vary based on the compiler used to build the B34S. On a Gateway P6 400 machine LF90 4.50i LF95 5.5b SOLVE time 9.718 9.22 DO time 41.69 13.73 Gain of SOLVE 4.3897 1.49 On a Dell Latitude 1MH machine using verions 8.10X the SOLVE and DO times were 2.6538162 and 3.8855877 respectively giving a gain of 1.4641510. LF95 7.10 was used to build b34s. LF90 appears to make a very slow DO loop!! LF95 is faster. In simple formulas the FORMULA and SOLVE commands are useful. With large complex sequences of commands, the DO loop cost may have to be "eaten" by the user since it is relative low in comparison to the cost of parsing the expression. A major advantage of the DO loop is that the logic is 100% clear is most cases. Speed can be increased by using variables for constants. This is because at parse time all scalars are made temps. Doing this outside the loop speeds things. Slow code do i=3,1000; x(i)=x(i)*2.; enddo; better code two=2.0; do i=3,1000; x(i)=x(i)*two; enddo; Vectorized code i=integers(3,1000); x=x(i)*2.; Compact vectorized code x=x(integers(3,1000))*2.; If all elements need to be changed the fastest code is x=x*2.; In the vectorized examples parse time is the same no matter whether there are 10 elements in x or 10,000,000. For speed gains from the use of masks, see # 18 below. Since B34S can create, compile and execute Fortran programs, for complex calculations a branch to Fortran is always an option. The larger the dataset the less the overhead costs. The Fortran example in matrix.mac illustrates dynamically building, compiling and calling a Fortran program from the matrix command. b34sexec matrix; call open(70,'_test.f'); call rewind(70); /$ 1234567890 call character(test," write(6,*)'This is a test # 2'" " n=1000 " " write(6,*)n " " do i=1,n " " write(6,*) sin(float(i)) " " enddo " " stop " " end "); call write(test,70); call close(70); /$ lf95 is Lahey Compiler /$ g77 is Linux Compiler /$ fortcl is script to run Lahey LF95 on Unix to link libs call dodos('lf95 _test.f'); * call dounix('g77 _test.f -o_test'); call dounix('lf95 _test.f -o_test'); * call dounix('fortcl _test.f -o_test'); call dodos('_test > testout':); call dounix('./_test > testout':); call open(71,'testout'); call character(test2, ' '); call read(test2,71); call print(test2); testd=0.0; n=0; call read(n,71); testd=array(n:); call read(testd,71); call print(testd); call close(71); call dodos('erase testout'); call dodos('erase _test.f'); call dounix('rm testout'); call dounix('rm _test.f'); b34srun; For a more advanced example, see FORTRAN_2 test case in matrix.mac. 16. DO, IF and DOWHILE nested statement limits Do loop, DOWHILE loop and IF statement termination should be hit. If this is not done, the max IF statement limit or DO statement limit can be exceeded depending on program logic. This "limitation" comes from having IF and DO loops outside programs, subroutines or functions. When the DO loop or IF structure starts to run the program does not know the end. The code loop continue; if(dabs(z1-z2).gt.1.d-13)then; z2=z1; z1=dlog(z1)+c; go to loop; endif; will never hit endif; The b34s parser will not know the position of this statement and the max IF statement limit could be hit if the IF structure was parsed many times. A better approach is not to use an IF structure in this situation. The correct code is: loop continue; if(dabs(z1-z2).le.1.d-13)go to nextstep; z2=z1; z1=dlog(z1)+c; go to loop; nextstep continue; 17. Mask Issues Assume an array x where for x < 0 we want y=x**2. while for x GE 0.0 we want y=2*x. A slow way to do this would be: do i=1,norows(x); if(x(i) .lt. 0.0)y(i)=x(i)**2.; if(x(i) .ge. 0.0)y(i)=x(i)*2. ; enddo; since the larger the X array the more parsing is required since the do loop cycles more times. A vectorized way to do the same calculation is to define two masks. Mask1 = 0.0 if the logical expression is false, = 1.0 if it is true. Fast code would be mask1= x .lt. 0.0 ; mask2= x .ge. 0.0 ; y= mask1*(x**2.0) + mask2*(x*2.0); Faster code would be y= (x .lt. 0.0)*(x**2.0) + (x .ge. 0.0 )*(x*2.0); Complete problem: b34sexec matrix; call print('If X GE 0.0 y=20*x. Else y=2*x':); x=rn(array(20:)); mask1= x .lt. 0.0 ; mask2= x .ge. 0.0 ; y= mask1*(x*2.0) + mask2*(x*20.); call tabulate(x,y,mask1,mask2); b34srun; Compact code (placing the logical calculation in expression) is: b34sexec matrix; call print('If X GE 0.0 y=20*x. Else y=2*x':); x=rn(array(20:)); y= (x.lt.0.0)*(x*2.0) + (x.ge.0.0)*(x*20.); call tabulate(x,y); b34srun; Logical mask expressions can be used in function and subroutine calls to speed calculation. 18. N Dimensional Objects While the Matrix command saves only 1 and 2 dimensional objects, it is possible to save and address n dimensional objects in 1 d arrays. n dimensional objects are saved by col. The commands: nn=index(i,j,k:); x=array(nn); call setndimv(index(i,j,k),index(1,2,3),x,value); will make an 3 dimensional(i,j,k) object x and place value in the 1 2 3 position. The function call yy=getndimv(index(i,j,k),index(1,2,3),x); or yy=x(index(i,j,k:1,2,3)); can be used to pull a value out. For example to define the 4 dimensional object x with dimensions 2 3 4 5 nn=index(2,3,4,5:); x=array(nn:); To fill this array with values 1.,...,norows(x) x=dfloat(integers(norows(x))); to set the 1 2 3 1 value to 100. call setndim(index(2,3,4,5),index(1,2,3,1),x,100.); Examples: b34sexec matrix; x=rn(array(index(4,4,4:):)); call print(x,getndimv(index(4,4,4),index(1,2,1),x)); do k=1,4; do i=1,4; do j=1,4; test=getndimv(index(4,4,4),index(i,j,k),x); call print(i,j,k,test); enddo; enddo; enddo; b34srun; b34sexec matrix; xx=index(1,2,3,4,5,4,3); call names(all); call print(xx); call print('Integer*4 Array ',index(1 2 3 4 5 4 3)); call print('# elements in 1 2 3 4 is 24',index(2 3 4:)); call print('Position of 1 2 in a 4 by 4 is 5',index(4 4:1 2):); call print('Integer*4 Array ',index(1,2,3,4,5 4 3)); call print('# elements in 1 2 3 5 is 30',index(2,3,5:)); call print('Position of 1 3 in a 4 by 4 is 9',index(4,4:1,3):); b34srun; b34sexec matrix; mm=index(4,5,6:); xx=rn(array(mm:)); idim =index(4,5,6); idim2=index(2,2,2); call setndimv(idim,idim2,xx,10.); vv= getndimv(idim,idim2 ,xx); call print(xx,vv); b34srun; N dimensional arrays are very good at saving data. What thgey are not good is directly processing the data. Here we need to extract what is needed into 1D or 2D objects, make the needed calculations and replace the values back in the structure. 19. Complex Math Issues If x=complex(1.5,1.5); y=complex(1.0,0.0); a=x*y; produces a=(1.5,1.5); to zero out imag part of a use a=complex(real(x*y),0.0); Assume x is a real*8 array and y=2., the statement newx=x**y; squares all elements. If y is not equal to a positive integer greater than 1. then x must be made complex prior to the calculation. Now assume x is a matrix. If y is a positive integer, then x can be real. If y is not an integer, then both x and y must be complex*16. The calculation proceeds as follows. First do a eigenvalue factorization of x as x = c * dlamda * inv(c) where c is the eigenvectors of x and dlamda is a diagonal matrix with eigenvalues along the diagonal. We note that x**y = c * (dlamda**y) * inv(c) which has transformed what would have been a very difficult calculation into an element by element operation along the diagonal. RG option: To speed up calculation x is inspected. If the imag part is all 0.0, the EISPACK RG routine was used. This is substantially faster than the EISPACK CG which is used when x is a complex number. The RG calculation seems to differ in accuracy from MATLAB. It appears that MATLAB is using a complex eigenvalue analyis to form the power. To make B34S run simular to MATLAB, using RG is no longer the default. If the user wants to force RG to be used in these situations where x is real place line b34sexec options debugsubs(b34smatpower2); b34srun; prior to calling b34sexec matrix; 20. Inversion Issues. Assume X is a n by n real*8, real*16 or complex*16 or complex*32 matrix or VPA matrix. The command invx=inv(x); will invert x using the LINPACK LU factorization routines. If x is not full rank, the progam will stop with an error message. The optional argument rcond in invx=inv(x,rcond); will return the condition. If it is known that x is symmetric or positive definate, the key words :smat or :pdmat can be used to speed up the calculation 3 or more times. The default is the LU factorization using LINPACK routine DGECO/DGEDI and ZGECO/ZGEDI. The optional argument :gmat will force use of the LAPACK routines DGETRF/ZGETRF and DGETRI/ZGETRI. The switch :pdmat2 will call the LAPACK routines DPOTRF/DPOTRI in place of LINPACK. LAPACK is not available for real*16, complex*32 or VPA matrices. For example invx=inv(x); invx=inv(x:gmat); invx=inv(x:smat); invx=inv(x:pdmat); invx=inv(x:pdmat2); A number of jobs can be run to illustrate the speed gains. default uses LINPACK DGECO/DGEDI and ZGECO/ZGEDI gmat uses LAPACK DGETRF/ZGETRF and DGETRI/ZGETRI smat uses LINPACK DSICO/DSIDI and ZSICO/ZSIDI pdmat uses LINPACK DPOCO/DPODI and ZPOCO/ZPODI pdmat2 uses LAPACK DPOTRF/DPOTRI and ZPOTRF/ZPOTRI The LAPACK library has routines DGESVX and ZGESVX that allows solutions of the system X*A=B where x is n by n, a is n by m and b is n by m using refinement and or balancing of the system. While it is not space efficient, if b is set as the identity matrix, such routines can be used to refine the inverse. The optional key words :refine and :refinee allow solution of the inverse where the solution is refined and or refined and equilibrated. invx=inv(x:refine); invx=inv(x:refinee); The LAPACK rcond estimator differs slightly from the LINPACK rcond estimator. For a symmetric positive definate matrix, the Cholesky routines r=pdfac(x); xinv=pdinv(r); can be used. Rank problems can be detected with optional arguments. r=pdfac(x,rcond); r=pdfac(x,rcond,ibad); and solutions of equations can be calculated with a=pdsolv(r,b); The Cholesky r can be updated and down dated using the commands pdfacdd and pdfacud. The commands pdfac, pdinv, pdsolv, pdfacdd and pdfacud work with real*4, real*8, real*16, complex*16 and complex*32 matrixes. While the above solvers all use LINPACK, it is well known that the QR method can be used to obtain r. Assuming x is n by k, the command r1 = qrfac(x); produces the same r as r2=pdfac(transpose(x)*x); Using the QR routines LINPACK DQRDC and ZQRDC routines, and their real*16 and complex*32 counterparts qqrdc and cqqrdc, it is possible to solve the OLS problem as: qr=qrfac(x,pivot); b=qrsolve(qr,pivot,y,info); where qrsolve uses the LINPACK DQRSL / ZQRSL or QQRSL and CQQRSL routines QRFAC and QRSOLV work for real*4,real*8, real*16, complex*16 or complex*32 data. Real and complex VPA (variable precision arimatic) For a real*8 matrix, the generalized inverse can be obtained by ginv=pinv(x,irank); ginv=pinv(x,irank,toll); where: irank = an estimate of the rank of x toll = tolerance that is used to set the singular values to zero. At present pinv works for a real*8 matrix x. IMSL routine DLSGRR is used for the calculation. In many applications one may want to determine easily if a matrix is full rank and take corrective action without placing the inverse in the code. The call call gmfac(x,l,u,info); Factors n by n matrix x such that x = L*U where L is lower triangular with 1.0 on the diagonal and U is upper triangular. In contrast to x=inv(xx); which uses LINPAC DGECO/DGEFA/DGEDI ZGECO/ZGEFA/ZGEDI, and for real*16 and complex*32 QGECO/QGEFA/QGEDI CQGECO/CQGEFA/CQGEDI GMFAC uses the LAPACK routines DGETRF and ZGETRF. GMFAC optionally will return info > 0 if U(i,i)=0. GMFAC will not run on real*16 and complex*32 data types. The commands call gminv(x,xinv); call gminv(x,xinv,info); call gminv(x,xinv,info,rcond); invert a general matrix using LAPACK. This code may be faster than inv( ) if the rank of the matrix is greater than 200 by 200. If the optional argument info is present, the routine will not stop if there is a problem. This allows the user to take automatic corrective action. The routine GMINV is the fastest large general matrix routine. It does not do refinement by default. The rcond option takes time and usually is not needed. The rcond value from LAPACK is not the same as LINPACK. The routine GMINV uses LAPACK DGETRF/ZGETRF and DGETRI/ZGETRI and optionally DGECON/DGECON if the condition is needed. The command call gmsolv(x,b,ans,info); where x is a n by n matrix, b a n by k matrix solves the system ans*x=b. If the optional argument info is present, the program will not stop if there is a rank problem. The LAPACK routines DGETRF/DGETRS and ZGETRF/ZGETRS are usually used. If the key words :refine or :refinee are present, the LAPACK routines DGESVX and ZGESVX will be used to refine the solution. This will take substantially more time. Unless the optional key words are used, gmsolv is the fastest way to solve the general system ans*x=b where x is greater than or equal to 200 by 200. The larger b, the more the gain in speed. The command s=svd(x,ibad,job,u,v); Calculates singular value decomposition of x which must be real*8 or real*16, complex*16 or complex*32. LINPACK routines DSVDC and ZSVDC are used for real*8 or complex*16 x. For real*16 and complex*32, these have been converted to QSVDC and CQSVDC respectively. For added detail on these options, please look at the help files and the example files. Inversion options under the matrix command have been designed to be both easy to use and flexible. If the series are real*16 or complex*32, only the LINPACK routines are available since LAPACK is not currently available for these data types. The inv( ) command also works for VPA data. Substantial accuracy gains can be obtained. For further detail, see sections 23 & 24 below. 21. Fortran and other Language Interfaces While the B34S Matrix language is "vectorized" and can easily handle most tasks, recursive systems that involve DO and DOWHILE loops are quite slow. An alternative is to code the desired calculation in the user language of choice and link to this program inside a B34S matrix command subroutine. While one possible implementation of this facility might be via a W2K DLL, this approach was rejected as being not portable and overly complex. The implementation discussed below allows the user to use any language to develop the program and communicate with the b34s processor with files. The below listed examples illustrate this technique. If there is sufficient user demand, at a later point a user DLL for the matrix command might be developed. The downside is that such an implementation would require the user directly access the matrix commmand named storage array. It is the view of the developer of B34S that only the most experienced programmers are capable of this complexity. A programmer without the proper training could easily "kill" or worse still, damage the B34S MATRIX command resulting in a hard-to-detect error being generated. The file IO approach will not allow direct access to the B34S matrix command arrays and is safer but slower. If most of the cpu time is in the calculation, not in passing the data, this cost will be minimal. Example of a user fortran program being called. The below listed example calls the program _test.exe on W2K which in turn writes a file of 1000 sin values. These are read into the matrix command. Except for the compile command this job is 100% portable between Linux and W2K. The user is NOT required to write in FORTRAN. Any user programming language or user external program is allowed provided it can read and make files for the IO into the B34S matrix command. Simple Example Program logic: 1. Open file _test.f 2. Copy lines in the 2-D character*1 array test into this file. 3. Compile the fortran. 4. Run the program _test and place output in testout 5. Open testout in the B34S matrix command and read the character first line. The number of elements line into valiable n. Allocate a 1d array testd of length n. Read data into testd. Print testd. b34sexec matrix; call open(70,'_test.f'); call rewind(70); /$ 1234567890 call character(test," write(6,*)'This is a test # 2'" " n=1000 " " write(6,*)n " " do i=1,n " " write(6,*) sin(float(i)) " " enddo " " stop " " end "); call write(test,70); call close(70); call dodos('lf95 _test.f'); call dounix('lf95 _test.f -o_test'); call dodos('_test > testout':); call dounix('./_test > testout':); call open(71,'testout'); call character(test2,' '); call read(test2,71); call print(test2); testd=0.0; n=0; call read(n,71); testd=array(n:); call read(testd,71); call print(testd); call close(71); call dodos('erase testout'); call dodos('erase _test.f'); call dounix('rm testout'); call dounix('rm _test.f'); b34srun; Complex Example Program Logic 1. Fortran program logic: Read data from data.dat into array data1. Read number of parameters in the model into npar. Read Model npar parameters into array parm. Calculate the function value and write to file testout. 2. Matrix program test writes the parameters as they change to a file and calls the compiled fortran program which will return with a new func value. Things to be aware of: The commands call dodos(' '); call dounit(' '); were not given in the form call dodos(' ':); call dounit(' ':); since there is no screen writting. By not using the : the flash is avoided. The matrix command write used for data was call write(y,72,'(3e25.16)'); in place of the more general call write(y,72); because the line got too long on Linux. The lines call out.... were added to show progress. These are only seen if the program is running in the Display Manager. b34sexec options ginclude('b34sdata.mac') member(lee4); b34srun; b34sexec matrix ; call loaddata ; * The data has been generated by GAUSS by with settings $ * a1 = GMA = 0.09 $ * b1_n = GAR = 0.5 ( When Negative) $ * b1 = GAR = 0.01 $ * call echooff ; /$ Setup fortran call open(70,'_test.f'); call rewind(70); call character(fortran, /$234567890 " implicit real*8(a-h,o-z) " " parameter(nn=10000) " " dimension data1(nn) " " dimension res1(nn) " " dimension res2(nn) " " dimension parm(100) " " call dcopy(nn,0.0d+00,0,data1,1)" " call dcopy(nn,0.0d+00,0,res2 ,1)" " open(unit=8,file='data.dat') " " open(unit=9,file='tdata.dat') " " read(8,*)nob " " read(8,*)(data1(ii),ii=1,nob) " " read(9,*)npar " " read(9,*)(parm(ii),ii=1,npar) " " read(9,*) res2(1) " " close(unit=9) " " " " do i=1,nob " " res1(i)=data1(i)-parm(3) " " enddo " " " " func=0.0d+00 " " do i=2,nob " " res2(i) =parm(1)+(parm(2)* res2(i-1) ) +" " * (parm(4)*(res1(i-1)**2) ) " " func=func+(dlog(dabs(res2(i))))+ " " * ((res1(i)**2)/res2(i)) " " enddo " " func=-.5d+00*func " " close(unit=8) " " open(unit=8,file='testout') " " write(8,fmt='(e25.16)')func " " close(unit=8) " " stop " " end "); call write(fortran,70); call close(70); maxlag=0 ; y=doo1 ; y=y-mean(y) ; * compile fortran and save data; call dodos('lf95 _test.f' ); * call dounix('lf95 _test.f -o_test'); call dounix('fortlc _test'); call open(72,'data.dat'); call rewind(72); call write(norows(y),72); call write(y,72,'(3e25.16)'); call close(72); v=variance(y) ; arch=array(norows(y):) + dsqrt(v); i=2; j=norows(y); * parm=array(:.0001 .0001 .0001 .0001); * parm(1)=v; * parm(3)=mean(y); count=0.0; call echooff; program test; call open(72,'tdata.dat'); call rewind(72); npar=4; call write(npar,72); call write(parm,72,'(e25.16)'); arch(1)=0.0d+00 ; call write(arch(1),72,'(e25.16)'); call close(72); call dodos('_test'); call dounix('./_test '); call open(71,'testout'); func=0.0; call read(func,71); call close(71); count=count+1.0; call outdouble(10,5 ,func); call outdouble(10,6 ,count); call outdouble(10,7, parm(1)); call outdouble(10,8, parm(2)); call outdouble(10,9, parm(3)); call outdouble(10,10,parm(4)); return; end; ll =array(4: -.1e+10, .1e-10,.1e-10,.1e-10); uu =array(4: .1e+10, .1e+10,.1e+10,.1e+10); rvec=array(4: .1 .1, .1, .1); parm=rvec; * call names(all); call cmaxf2(func :name test :parms parm :ivalue rvec :maxit 2000 :maxfun 2000 :maxg 2000 :lower ll :upper uu :print); *call dodos('erase testout'); *call dodos('erase _test.exe'); *call dounix('rm testout'); *call dounix('rm _test'); b34srun; 22. Polynomial Matrix Operations A polynomial matrix is one where the coeffiients are themselves coefficients. A two by two system with max order 3 could be saved in form B(1,1)(L) B(1,2)(L) B(2,1)(L) B(2,2)(L) where the zero order terms (constant) are saved in the term. The VAREST command saves the matrix in a form where the constant is placed last. This is NOT :BYVAR format which would have occured if call olsq(y y{1 to maxlag} x{1 to maxlag} were used. B(1,1)(L) B(1,2)(L) c(1) B(2,1)(L) B(2,2)(L) c(2) For a two variable system with 6 lags the coefficients are saved in a 2 by ((2*6)+1) matrix. VAREST also saves [I - B(L)] which when inverted gives the phi weights. The matrix polynomial routines POLYMINV (used to calculate the inverse of a matrix), POLYMMULT (used to calculate a multiplication) and POLYMDISP (used to print or display such a matrix) use the convention: B(1,1)(0) B(1,2)(0) B(1,1)(k) B(1,2)(k) B(2,1)(0) B(2,2)(0) B(2,1)(k) B(2,2)(k) which is termed :BYORDER The command POLYMCONV can be used to convert from one system to another. Each matrix in :BYORDER form has a three element integer variable of the form (nrow,ncol,degree+1) that allows it to be decoded. If nrow=ncol then the matrix can be inverted. A matrix in :BYVAR form does not need an index to decode it but must be saved as a n by ((n*order)+1) system. For example a 2 variable VAR with mar order 5 could be saved as a 2 by 11. The first row would be 5 lags of x(1) and 5 lags x(2) plus the constant. For further detail see the VAREST command example. A simple setup is: b34sexec matrix; call loaddata; call load(buildlag); call load(varest); call echooff; ibegin=1; iend=296; nlag=2; nterms=10; x=catcol(gasin,gasout); call varest(x,nlag,ibegin,iend,beta,t,sigma,corr,residual,1, a,ai,varx,varxhat,rsq); call print(beta,t,sigma,corr); call tabulate(varx,varxhat,rsq); call polymdisp(:display a ai); call polyminv(a ai psi ipsi nterms); call polymdisp(:display psi ipsi); b34srun; 23. Real*16 and Complex*32 Data types While the usual data types in the MATRIX command are real*8 and complex*16, there is a limited real*16 and complex*32 capability. In addition there are improvements to accuracy in real*8, real*16, complex*16 and complex*32 by modification the the BLAS routines. If higher accuracy is needed, the VPA option, discussed in the next section, can be used. Examples: Option 1. Works for all commands that use routines that involve BLAS routines. The command call real16on(:real16math); will perform math with real*8 objects using real*16 math and complex*16 objects using complex*32 math. The command call real16on; uses the IMSL routines dqadd and dqmult to increase accuracy. It is faster than call real16on(:real16math). These options can be turned off with the command call real16off; For example the default test case takes 2.04, with real16on it takes 3.14 while, with real16on(:real16math) set, the time is 5.95. In most cases this added accuracy is not needed. For real*16 and complex*32 objects, the command call real32on; uses added accuracy for qdot, qsum, qasum, qaxpy and qscal and added accuracy for cqdot, cqsum, cqasom, cqaxpy and cqscal call real32off; turns off this feature. The command call real32_vpa; uses vpa math to increase real*16 accuracy. This will run slowly but variable precision math is supported. This allows vpa math to support a number of otherwise real*16 commands. At a later date this may enhance complex*32 objects. For now it is experimental. Option 2. Real*16 and complex*32 variables can be created from real*8 and complex*16 variables with the r8tor16 and c16toc32 commands. The commands r16tor8 and c32to16 can be used to convert series back to the default datatypes. A series can be created to be a specific data type with the command kindas( ). Assume x is real*8 and xr16 is real*16. The user wants one_r8 to be real*8 and oner_16 to be real*16. The following commands can be used. one_r8=kindas(x,1.0); oner_16=kindas(xr16,1.0); A more accurate way to proceed is oner_16=kindas(xr16,real16('1.0')); oner_r8=kindas(x, real16('1.0')); While array and matrix math is 100% available, only a subset of commands are available for real*16 and complex*32. More detail is available under the help commands of the specific commands. See the matrix example r16c32 and the examples math5 and math6 for applications. As of 24 June 2003, INV, EIG, SEIG, SVD and a number of production commands such as DEXP, DLOG, SUMSQ, SUM etc have been converted. While LINPACK, and EISPACK have been converted to support real*16 and complex*32, FFTPACK and LAPACK have not. Warning: While real*8 and compolex*16 data can be converted to real*16 and complex*32 with the commands r8tor16( ) and c16toc32 respectively, much accuracy will be lost from the initial data read into real*8 and complex*16. For more accuracy it is it necessary to read into real*16 or complex*32 directly. For an example, of how this is done see the Filippelli (filippelli.dat) data which is in stattest.mac. 24. Variable Precision Arithmetic. The B34S has the capability of variable precision math where the number of digits calculated is up to 1780. For more detail on this facility see the help documents for subroutine VPASET and function VPA that provide the interface into this facility. Full matrix and array math is provided as well as the inverse or real and complex VPA objects. Unlike a number of other software systems, this facility is 100% integrated in the b34s system. For example, once a vpa number is built with the vpa(function), calculations can be made in the usual manner as will be illustrated below. Real, integer and complex objects have been implemented. The facility allows hugh numbers, both real and integer, to be used in calculations. Given the default settings of MBASE is 10**7 and NDIGMX = 256, integers less than 10**1792 can be used. In contrast, the range for integer*4 is -2,147,483,648 - 2,147,483,647. For integer*8 these become -9,223,372,036,854,775,808 - 9,223,372,036,854,775,807. This facility was built with the FM_ZM Library version 1.2 built by David M. Smith. The basic reference is: Algorithm 693, ACM Transactions on Mathematical Software, Vol. 17, No. 2, June 1991, pages 273-283. although a newer version of the libraty was obtained from the web. The below listed program calculates 2.0/4.11 a number of ways and shows the various accuracy values obtained. b34sexec matrix; call echooff; * Accuracy differences depending on data precision; call print(' 2.0/4.11 using different precisions':); call fprint(:clear :col 32 :string ' 10 20 30 40 50' :print :clear :col 32 :string '12345678901234567890123456789012345678901234567890' :print); call fprint(:clear :col 1 :string 'str => vpa' :col 30 :display vpa('2. ')/vpa('4.11') 'e70.52' :print :clear :col 30 :display vpa('2.00')/vpa('4.11') 'e50.32' :print); call fprint(:clear :col 1 :string 'real*8 => vpa' :col 30 :display vpa(2.00)/vpa(4.11) 'e50.32' :print); call fprint(:clear :col 1 :string 'real*8 => real*16' :col 18 :display r8tor16(2.)/r8tor16(4.11) '(e50.32)' :print); call fprint(:clear :col 1 :string 'str => real*16' :col 18 :display real16('2.00')/real16('4.11') '(e50.32)' :print); call fprint(:clear :col 1 :string 'real*8 => real*8' :col 18 :display array(:2.00/4.11) '(e50.32)' :print); call fprint(:clear :col 1 :string 'real*4 => real*4' :col 18 :display sngl(2.00)/sngl(4.11) '(e50.32)' :print); b34srun; The Above job shows 2./4.11 using various precisions. The VPA result to 52 digits is: .4866180048661800486618004866180048661800486618004866M+0 and shows the repeating pattern that will continue. When the data (2.0 and 4.11) was read with a string input, no accuracy was lost with real*16 given the format e50.32. This degree of accuracy was not obtained with real*16 if the data was read as real*8 first and then converted to real*16. The above job was edited to display only 72 columns and is listed below. 2.0/4.11 using different precisions 10 20 30 40 1234567890123456789012345678901234567890 str => vpa .4866180048661800486618004866180048661800 .48661800486618004866180048661800M+0 real*8 => vpa .48661800486618001080454992664677M+0 real*8 => real*16 0.48661800486618001080454992664677E+00 str => real*16 0.48661800486618004866180048661800E+00 real*8 => real*8 0.48661800486618000000000000000000E+00 real*4 => real*4 0.48661798200000000000000000000000E+00 Even for this simple case, the real*4 result began to differ at the 6th digit (although with rounding the failure was at 10). The real*8 result failed at digit 17 as expected. When the real*16 data input was by string there was no loss of accuracy up to the 32nd digit. However, when the data was first read as real*8 and converted to real*16 accuracy problems emerged by the 17th digit. The calculation was limited by the precision of the data read. This finding was verified when we moved real*8 data into VPA and got the same result as the real*8 into real*16. The str into vpa is the "true" answer. The next example show inverse gains: /; /; Shows gains in accuracy of the inverse with vpa /; b34sexec matrix; call echooff; n=3; call vpaset(:ndigits 1750); x=rn(matrix(n,n:)); r16x=r8tor16(x); vpax=vpa(x); call print('Real*4 tests',sngl(x),inv(sngl(x)),sngl(x)*inv(sngl(x))); call print('Real*8 tests',x, inv(x), x*inv(x)); call print('Real*16 tests',r16x,inv(r16x),r16x*inv(r16x)); call print('VPA tests',vpax,inv(vpax),vpax*inv(vpax)); b34srun; Edited output produces: Real*4 tests Matrix of 3 by 3 elements (real*4) 1 2 3 1 2.05157 1.27773 -1.32010 2 1.08325 -1.22596 -1.52445 3 0.825589E-01 0.338526 -0.459242 Matrix of 3 by 3 elements (real*4) 1 2 3 1 0.521061 0.675528E-01 -1.72205 2 0.179445 -0.402323 0.819689 3 0.225948 -0.284424 -1.88285 Matrix of 3 by 3 elements (real*4) 1 2 3 1 1.00000 0.919681E-08 0.759197E-07 2 0.585772E-07 1.00000 0.822609E-07 3 0.878866E-08 0.560976E-08 1.00000 Real*8 tests X = Matrix of 3 by 3 elements 1 2 3 1 2.05157 1.27773 -1.32010 2 1.08325 -1.22596 -1.52445 3 0.825589E-01 0.338525 -0.459242 Matrix of 3 by 3 elements 1 2 3 1 0.521061 0.675528E-01 -1.72204 2 0.179445 -0.402323 0.819689 3 0.225948 -0.284424 -1.88285 Matrix of 3 by 3 elements 1 2 3 1 1.00000 -0.555112E-16 0.00000 2 0.00000 1.00000 0.00000 3 -0.138778E-16 -0.277556E-16 1.00000 Real*16 tests R16X = Matrix of 3 by 3 elements (real*16) 1 2 3 1 2.05157 1.27773 -1.32010 2 1.08325 -1.22596 -1.52445 3 0.825589E-01 0.338525 -0.459242 Matrix of 3 by 3 elements (real*16) 1 2 3 1 0.521061 0.675528E-01 -1.72204 2 0.179445 -0.402323 0.819689 3 0.225948 -0.284424 -1.88285 Matrix of 3 by 3 elements (real*16) 1 2 3 1 1.00000 -0.144445E-33 0.00000 2 0.00000 1.00000 0.00000 3 0.00000 -0.240741E-34 1.00000 VPA tests VPAX = Matrix of 3 by 3 elements VPA - FM 1 2 3 1 .205157M+1 .127773M+1 -.132010M+1 2 .108325M+1 -.122596M+1 -.152445M+1 3 .825589M-1 .338525M+0 -.459242M+0 Matrix of 3 by 3 elements VPA - FM 1 2 3 1 .521061M+0 .675528M-1 -.172205M+1 2 .179445M+0 -.402323M+0 .819689M+0 3 .225948M+0 -.284424M+0 -.188285M+1 Matrix of 3 by 3 elements VPA - FM 1 2 3 1 .100000M+1 -.139488M-1757 .000000M 0 2 .168219M-1758 .100000M+1 .000000M 0 3 -.138187M-1758 .173950M-1758 .100000M+1 and shows the worse accuracy of the off diagonal of the real*4 case was 0.822609E-07 while the worse VPA error was .168219e-1758. The problem consists of random numbers and is designed to be easy. There can be gains is the data is read intro character then into higher precision. the job VPA8 in matrix,mac shows this and provides a template for character data reading: VPA8 Very Hard OLS Data. /; /; Problem very very hard!!! /; Illustrates reading into real*8 and then into real*16 and VPA /; Reading from character into real*16 /; Reading from chapter into VPA /; b34sexec data heading('Bruce McCullough Test Case'); input y x; datacards; 10000000001 1000000000.000 10000000002 1000000000.000 10000000003 1000000000.900 10000000004 1000000001.100 10000000005 1000000001.010 10000000006 1000000000.990 10000000007 1000000001.100 10000000008 1000000000.999 10000000009 1000000000.000 10000000010 1000000000.001 b34sreturn; b34srun; b34sexec matrix; call loaddata; call echooff; /; /; Data reread into character to see what happens /; datacards; 10000000001 1000000000.000 10000000002 1000000000.000 10000000003 1000000000.900 10000000004 1000000001.100 10000000005 1000000001.010 10000000006 1000000000.990 10000000007 1000000001.100 10000000008 1000000000.999 10000000009 1000000000.000 10000000010 1000000000.001 b34sreturn; y16=r8tor16(y); x16=r8tor16(x); /; Cholesky on real*8 fails /; call olsq(y x :print); /; QR fails /; call olsq(y x :qr :print); /; real*16 Cholesky fails /; call olsq(y16 x16 :print); call olsq(y16 x16 :qr :print :savex); coef8_16=rollright(%coef); /; Go to vpa 64 digits and 20 digits to see if makes a difference y_vpa=vpa(%y); x_vpa=vpa(%x); beta_vpa=inv(transpose(x_vpa)*x_vpa)*transpose(x_vpa)*y_vpa; call print('Beta from VPA - default 64 digits accuracy ',beta_vpa); call print(beta_vpa(1)); call print(beta_vpa(2)); call vpaset(:ndigits 200); y_vpa=vpa(%y); x_vpa=vpa(%x); beta_vpa=inv(transpose(x_vpa)*x_vpa)*transpose(x_vpa)*y_vpa; call print('Beta from VPA - 200 digits accuracy ',beta_vpa); call print(beta_vpa(1)); call print(beta_vpa(2)); /; Reading data as character into VPA and real*16 /; Makes a difference !!!!!!!!!!!!!!!!!!!!!!!!!!! data16=r8tor16(array(20:)); call read(data16,4); call load(ntokin :staging); call load(getr16 :staging); call load(getvpa :staging); call load(getchar :staging); call vpaset(:ndigits 200); vpdata=vpa(data16); vpx=vpdata*vpa('0.0'); call rewind(4); cc=c1array(10000:); call read(cc,4); call ntokin(cc,nfind,0,ibad); call print(nfind ); call getvpa(cc,nfind,vpx,ibad); call free(cc); xm=matrix(10,2:data16); y16=array(:xm(,1)); x16=matrix(10,2:data16); x16(,1)=real16('1.'); /; This will not solve due to condition /; call olsq(y16 x16 :print :noint); /; call print('Here data read into real*16 directly - Note difference':); call olsq(y16 x16 :qr :print :noint); call print(' ':); call print('Display More digits to see what real*16 has':); call print(' ':); call fprint(:clear :display %coef(1) '(1pe48.32)' :print); call fprint(:clear :display %coef(2) '(1pe48.32)' :print); call print(' ':); call print('Compare with reading into real*8 then to real*16':); call print(' ':); call fprint(:clear :display coef8_16(1) '(1pe48.32)' :print); call fprint(:clear :display coef8_16(2) '(1pe48.32)' :print); call print(' ':); /; now run with vpa vpxm=matrix(10,2:vpx); vpy= vector(:vpxm(,1)); vpx=vpxm; vpx(,1)=vpa('1.0'); betavpa=inv(transpose(vpx)*vpx)*transpose(vpx)*vpy; call print(' ':); call print('Beta from VPA where character read into VPA':); call print('This is the right way to attack the problem':); call print('For gain compare with real*16 where character read':); call print(betavpa(1)); call print(betavpa(2)); call print(' ':); call print('Beta from VPA where vpa converted from real*8':); beta_vpa=rollright(beta_vpa);; call print(beta_vpa(1)); call print(beta_vpa(2)); call print(' ':); call print('Difference of coef # 1 and # 2':); call print('Shows effect of not reading in correct precision':); diffcoef=vpa(vector(2:)); diffcoef(1)=betavpa(1)-beta_vpa(1); diffcoef(2)=betavpa(2)-beta_vpa(2); call print(diffcoef(1)); call print(diffcoef(2)); b34srun; Heavily Edited output is: Variable # Cases Mean Std Deviation Variance Y 1 10 0.1000000001E+11 3.027650354 9.166666667 X 2 10 1000000001. 0.5278048818 0.2785779932 CONSTANT 3 10 1.000000000 0.000000000 0.000000000 Number of observations in data file 10 Current missing variable code 1.000000000000000E+31 B34S(r) Matrix Command. d/m/y 29/ 7/07. h:m:s 8:37:30. Ordinary Least Squares Estimation - Real*16 & QR Variable Lag Coefficient SE t X16 0 0.96522007E-01 2.0278033 0.47599294E-01 CONSTANT 0 0.99034780E+10 0.20278033E+10 4.8838454 Beta from VPA - default 64 digits accuracy BETA_VPA= Vector of 2 elements VPA - FM .965220M-1 .990348M+10 9.65220067890258342009546806150744940551143108123572854800000M-2 9.90347799865209574142761316341563647882994345591008434102590M+9 Beta from VPA - 200 digits accuracy BETA_VPA= Vector of 2 elements VPA - FM .965220M-1 .990348M+10 9.6522006789025834200954680615074494055114312994576814995573154883862781 971061687532225195211872707883443648855216365597479743619588764846877171 146998987325841273459533669904135897504647183082754700000M-2 9.9034779986520957414276131634156364788299434544608648119274223099391811 078965853733763483342081454150397522719740321509191330273731665955348598 768449388118628589173691417115419955283479769263019114930M+9 Display More digits to see what real*16 has 9.90347806584471095667600829609938E+09 9.65219395964106601701812011968070E-02 Compare with reading into real*8 then to real*16 9.90347799865209574142752881434979E+09 9.65220067890258342010390296808680E-02 Beta from VPA where character read into VPA This is the right way to attack the problem For gain compare with real*16 where character read 9.9034780658447109566760077568540548388203264036962318951564333468144967 976254007455322706347553966533211125389976555538803813972707424451639716 305267784566221628731949001317005969203917343089228659810M+9 9.6521939596410660170181740442134299509971673602685383945928568978486775 297722321536118749107570909723269206071150230416216962175365207909055592 648697631862131571369199609764191317652106807309346900000M-2 Beta from VPA where vpa converted from real*8 9.9034779986520957414276131634156364788299434544608648119274223099391811 078965853733763483342081454150397522719740321509191330273731665955348598 768449388118628589173691417115419955283479769263019114930M+9 9.6522006789025834200954680615074494055114312994576814995573154883862781 971061687532225195211872707883443648855216365597479743619588764846877171 146998987325841273459533669904135897504647183082754700000M-2 Difference of coef # 1 and # 2 Shows effect of not reading in correct precision 6.7192615215248394593438418359990382949235367083229011036875315689728815 372155922300547251238281360267023623402961248369897575849629111753681839 644759303955825758420158601392043757382620954488000000000M+1 -6.719261517403077294017294019454514263939189143104964458590537600667333 936599610644610430179816017444278406613518126278144422355693782157849830 1355463709702090334060139944579852540375773407800000000000M-8 This example shows that reading into real*8 first then converting to real*16 and or VPA causes a loss of accuracy in difficult problems. 25. Order of Operation issues. The command call bestreg (y x1 x2 x3 x4 :crit -3 :print); will not work as intended since the parser will attempt the operation crit-3. There are two ways to proceed. One it to pass a variable mthree=-3; call bestreg (y x1 x2 x3 x4 :crit mthree :print); A more elegant way is to force creation of a temp variable call bestreg (y x1 x2 x3 x4 :crit sfam(-3) :print); 26. Power of a matrix, array and a vector If x is an array the statements x=array(2:1 2 3); xnew=x**2.; produce an array of 1 4 9 If x is defined as a vector, x=vfam(x); xvnew =x**2.; produces 14 or the same as x*x If x is a 2-D array, then x=array(3,3:1 2 3 4 5 6 7 8 9) xnew=x*2 gives a matrix 1 4 9 16 25 36 49 64 81 If x is a matrix, then x_matnew = x**3.; is the same as x*x*x However the statement x_matnew=x**3.5; cannot be done this way. The solution relies on decomposing x as v*s*inv(v) where v are the eigenvectors of x and s is a diagonal matrix with the eigenvalues (e) of x along the diagonal. If x is a general matrix (not positive definate) these are possibly complex. The calculation for x_matnew is x_matnew = v * ss * inv(v) where ss is a diagonal matrix where ss(i,i)=e(i)**3.5. Some of these ideas are illustrated in the example file math_power in matrix.mac which is shown below. /; /; Matrix, array and Vector Power tests /; b34sexec matrix; aa=array(4:1 2 3 4); va=vfam(aa); aa16=r8tor16(aa); va16=r8tor16(va); caa=complex(aa); cva=complex(va); qcaa=qcomplex(aa16); qcva=qcomplex(va16); saa=sngl(aa); sva=vfam(saa); vpa_aa=vpa(aa); vpa_va=vpa(va); saa=sngl(aa); sva=vfam(saa); /; /; array and vector tests /; call print(saa, saa*saa, saa**sngl(2.)); call print(sva, sva*sva, sva**sngl(2.)); call print(aa16, aa16*aa16, aa16**r8tor16(2.)); call print(va16, va16*va16, va16**r8tor16(2.)); call print(caa, caa*caa, caa**complex(2.)); call print(cva, cva*cva, cva**complex(2.)); call print(qcaa, qcaa*qcaa, qcaa**qcomplex(r8tor16(2.)), qcva, qcva*qcva, qcva**qcomplex(r8tor16(2.))); call print(vpa_aa ,vpa_aa *vpa_aa, vpa_aa**vpa(2.), vpa_va ,vpa_va *vpa_va, vpa_va**vpa(2.)); call print('2 D Tests ++++++++++++++':); aa=array(3,3:1. 2. 30. 4. 5. 60. 7. 8. 90.); va=vfam(aa); aa16=r8tor16(aa); va16=r8tor16(va); caa=complex(aa); cva=complex(va); qcaa=qcomplex(aa16); qcva=qcomplex(va16); saa=sngl(aa); sva=vfam(saa); vpa_aa=vpa(aa); vpa_va=vpa(va); saa=sngl(aa); sva=vfam(saa); /; /; array and matrix /; call print(saa, saa*saa , saa**sngl(2.)); call print(sva, sva*sva , sva**sngl(2.)); call print( aa, aa* aa , aa**2.); call print( va, va* va , va**2.); /; /; This calculation will draw an error /; /; call print( va**2.5); /; call print(aa16, aa16*aa16, aa16**r8tor16(2.)); call print(va16, va16*va16, va16**r8tor16(2.)); call print(caa ,caa *caa, caa**complex(2.)); call print(cva ,cva *cva, cva**complex(2.)); call print(cva ,cva *cva*cva cva**complex(3.)); call print(qcaa, qcaa*qcaa, qcaa**qcomplex(r8tor16(2.))); call print(qcva, qcva*qcva, qcva**qcomplex(r8tor16(2.))); call print(vpa_aa ,vpa_aa *vpa_aa, vpa_aa**vpa(2.)); call print(vpa_aa ,vpa_aa *vpa_aa*vpa_aa, vpa_aa**vpa(3.)); call print(vpa_va ,vpa_va *vpa_va, vpa_va**vpa(2.)); call print(vpa_va ,vpa_va *vpa_va*vpa_va, vpa_va**vpa(3.)); /; /; this is an error since no eig capabilkity for vpa /; call print(vpa_va ,vpa_va *vpa_va*vpa_va, vpa_va**vpa(3.1)); /; /; Non integer power for complex*16 & complex*32 /; call print(eig(cva),eig(qcva)); call print('Non integer power complex*16', cva**complex(3.5)); call print('Non integer power complex*32', qcva**qcomplex(r8tor16(3.5))); /; /; Slow way to get an inverse /; * Test case for Real Matrix from IMSL Math (10) pp 295-297; * using eig( ) analysis to get inverse; x=complex(matrix(3,3:8.,-1.,-5.,-4., 4.,-2.,18.,-5.,-7.)); call print('X Matrix',x); call print(inv(x), x**complex(-1.)); call print(inv(c16toc32(x)), c16toc32(x)**c16toc32(complex(-1.))); b34srun; 26. Level issues The default command level is 100. If a subroutine is called the level is set to 101. If another routine is called, the level becomes 102. LEVEL=1 is the global level. Variables at the global level must be moved to the local level to be used because of assingment level issues. For example xlocal=xglobal; makes a local copy of xglobal. As an alternative, the command call makeglobal(x) Puts X at level 1. X can be passed to a subroutine where a local copy is automatically made. Any changes are placed back at the global level. As another alternative, the command call makelocal(x); will make x local. As notes, if X is global, statements such as y=x*x; are not allowed since it is not clear at what level to place the answer. If x is desired to be used, make a local copy. /; /; Global variable moved to a subroutine and made local. /; Calculations made and the result put back at global level. /; b34sexec matrix; y=array(:10. 20.); call makeglobal(y); call print(level(v)); call print(level(y)); subroutine useit(x); call names(all); call print(x); x=x+100.; return; end; call useit(y); call names(all); call print(y); b34srun; Matrix Programming Language key words CALL Call a subroutine call print('This is a string'); calls the print command. B34S MATRIX commands that do not return an argument across an equals are executed by the CALL sentence. The CALL sentence first looks in named storage for a routine with this name. If this is not found, then the built-in routines are used. While it is possible to have a user routine with the same name as a built in routine, this is not a good idea and can caluse unexpected errors. CONTINUE go to statement name continue; is a go to statement target. Note name must be le 8 characters Example: if(x.eq.0.0)go to test; call print('x is greater than 0.0'); test continue; DO Starts a do loop do i=1,10; Begins a DO loop. An alternative is: do i=1,10,1; Another alternative is for i=1,10; Notes on do loops: The code: do i=1,10; do j=i,10; s(i,j)=b(i)*a(j); enddo; enddo; is valid. As written, both do loops will completely execute. If in place of the above, the code had been: do i=1,10; do j=i,10; s(i,j)=b(i)*a(j); if(s(i,j).le..0001)go to done; enddo; done continue; call print('This is the end of loop one'); enddo; the results would not be as intended since B34S does not know to return to the statement do i=1,10; and will instead try to continue with the statement do j=i,10; since it may not "know" the exact location of the inner loop enddo; statement because it may not have found it. The solution is to manually terminate the inner loop by setting j outside the range. Corrected code would be: do i=1,10; do j=i,10; s(i,j)=b(i)*a(j); if(s(i,j).le..0001)go to done; enddo; done continue; j=0; call print('This is the end of loop one'); enddo; Example: do i=1,n; x(i)=y(i+1); enddo; A faster code would be i=integers(1,n); j=i+1; x(i)=y(j); Notes on DO loop variables. In Fortran it is not recommended that DO loop variables be used outside the loop. In B34S b34sexec matrix; do i=1,2; call print('in loop i = ',i); enddo; call print('out of loop i =,i); b34srun; Produces i=3 at the "out of loop" position. This is the same as Fortran. To check compile c test fortran do i=1,2 write(6,*)'In loop i was ',i enddo write(6,*)'Out of loop i was ',i end Note how closely B34S follows the Fortran standard and the language. DOWHILE Starts a dowhile loop dowhile(x.gt.0.0); Starts a DOWHILE loop. Example: b34sexec matrix; sum=0.0; add=1.; count=1.; tol=.1e-6; dowhile (add.gt.tol); oldsum=sum; sum=oldsum+((1./count)**3.); count=count+1.; add=sum-oldsum; enddowhile; call print('Sum was ',sum:); call print('Count was ',count); b34srun; Warning: Be sure you do not have an infinate loop. ENDDO ENDS a DO loop enddo; Ends a DO loop. next i; Can be used in place of enddo; Example: do i=1,n; x(i)=y(i+1); enddo; ENDDOWHILE ENDS a dowhile loop enddowhile; Ends a DOWHILE loop. Example: b34sexec matrix; sum=0.0; add=1.; ccount=1.; count=1.; tol=.1e-6; /$ outer dowhile does things 2 times call outstring(2,2,'We sum until we can add nothing!!'); call outstring(2,4,'Tol set as '); call outdouble(20,4,tol); call echooff; dowhile(ccount.ge.1..and.ccount.le.3.); sum=0.0; add=1.; count=1.; dowhile(add.gt.tol); oldsum=sum; sum=oldsum+((1./count)**3.); count=count+1.; call outdouble(2,6,add); add=sum-oldsum; enddowhile; ccount=ccount+1.; call print('Sum was ',sum:); call print('Count was ',count); call compress(1000); enddowhile; b34srun; END End of a program, function or Subroutine. end; Ends a PROGRAM, SUBROUTINE or FUNCTION. Example: subroutine test(x); call print('The mean of x is',mean(x)); return; end; EXITDO Exit a DO loop exitdo; Exits a do loop. Example: b34sexec matrix; call echooff; do j=1,4; do i=1,10; if(i.eq.8)exitdo; if(i.ge.0)then; call print('at 1 in if i was ',i:); if(i.ge.4)exitif; call print('at 2 in if Should never be gt 3 i was ',i:); endif; call print('in do loop ',i); enddo; enddo; b34srun; EXITIF Exit an IF loop exitdo; Exits a do loop. Example: b34sexec matrix; call echooff; do j=1,4; do i=1,10; if(i.eq.8)exitdo; if(i.ge.0)then; call print('at 1 in if i was ',i:); if(i.ge.4)exitif; call print('at 2 in if Should never be gt 3 i was ',i:); endif; call print('in do loop ',i); enddo; enddo; b34srun; FORMULA Defive a recursive formula. formula x=y(t-1); Defines a formula. Formulas are only executed for observation t in a SOLVE statement. Formula definations are saved at level 2. Formula results are saved at the current level and can be used in the usual analytical statements provided that they have been executed by a solve statement. More extensive help in given later in the help document. Brief example: * archvar and resid start out as variables =0.0; * formula statement updates ONLY for obs t; * For t=2, this value of u is used to get archvar; archvar=array(norows(y):); resid =array(norows(y):); formula archvar = a0 + a1 * (resid(t-1)**2.) ; formula resid=y(t) - b1 - b2*x1(t) - b3*dsqrt(archvar(t)); solve( archlogl=(-.5)*(dlog(archvar(t))+((resid(t)**2.)/archvar(t)) :range 2 norows(y) :block archvar resid); FOR Start a DO loop for i=1,10; Alternate do loop setup. Example: for i=1,n; x(i)=y(i+1); next i; is the same as do i=1,n; x(i)=y(i+1); enddo; NEXT i End of do loop next i; Alternate end of DO loop Example: for i=1,n; x(i)=y(i+1); next i; GO TO Transfer statement go to n; Transfers control to statement n CONTINUE; Note n must be le 8 characters Example: if(x.eq.0.0)go to test; call print('x is greater than 0.0'); test continue; Note: Statements such as if(k.eq.0)then; * statements here ; if(jj.gt.0)go to done; endif; should be avoided since the statement endif; will not be found and the # of if statements will be exceeded. Better code is bad=0; if(k.eq.0)then; * statements here; bad=1; endif; if(bad.eq.1)go to done; or if(k.eq.0)then; * statements here; if(jj.gt.0)exitif; endif; FUNCTION Beginning of a function. function somename(args); is the first line of a user function. Functions can have functions as arguments and themselves can be used as arguments. The command call compress; will be ignored if found in a function or a subroutine or a program called by a running function. Note: Repeated calls to user functions use up parse space. Thus for heavy computing, user subroutines should be employed in place of user functions. If functions are used outside subroutines, the option adjust_parse can be used. This "undocumented" command will NOT work if a function is used as an argument inside a subroutine call. Hence the use of the option adjust_parse should not be generally employed. As a second place alternative, the parse tree can be expanded by use of b34sexec options maxnsent(600) maxntoken(5000); b34srun; Examples: function tt(y); t=sum(y); return(t); end; testmean=tt(y)/dfloat(norows(y)); call mysub(tt(y),x,z); IF( ) Beginning of an IF structure Note: The WHERE statement operates on non scalar objects while the IF statement operates on scalar objects. if(x.eq.9)y=dabs(x); Simple IF statement. X must be a scalar. If the mask JJ (must be 0.0 or 1.0) was set the code: jj(1)=0.0; jj(2)=1.0; if(jj(i))call print('i was 2'); can be used. The commands NLPMIN1, NLPMIN2 and NLPMIN3 use mask technology. The statements: if(x.eq.9)then; call stop; endif; will STOP the program and get out of the MATRIX command if x=9. Note that IF statements ( ) must resolve to 0.0 or 1.0. A check for 0.0 of the form if(x.ne.0.0)y=p/x; will work in versions of B34S since November 7, 2004. An alternative is: if(x.ne.0.0)then; y=p/x; endif; MASKS are a feature related to IF statements but much faster. Since a logical expression resolves to be 0.0 or 1.0, a mask can be built with a expression such as mask = x .gt. 1.0 ; Masks are a way to vectorize what would be an IF. For tyeh below listed example note that x is an array, not a vector. For example: b34sexec matrix; call print('If X GE 0.0 y=20*x. Else y=2*x':); x=rn(array(20:)); y= (x.lt.0.0)*(x*2.0) + (x.ge.0.0)*(x*20.); call tabulate(x,y); b34srun; will run faster than b34sexec matrix; call print('If X GE 0.0 y=20*x. Else y=2*x':); x=rn(array(20:)); do i=1,norows(x); if(x(i).lt.0.0)y(i)=x(i)*2.0; if(x(i).ge.0.0)y(i)=x(i)*20.; enddo; call tabulate(x,y); b34srun; due to DO loop overhead. Note: Statements such as if(k.eq.0)then; * statements here ; go to done; endif; should be avoided since the statement endif; will not be found and the # of if statements will be exceeded. Better code is bad=0; if(k.eq.0)then; * statements here; bad=1; endif; if(bad.eq.1)go to done; Examples of code to change an element of an array using a mask. x=integers(1,10); xx=x; /$ /$ Note that Where sets one element to 99 rest to 0.0; /$ where(x.eq.5)x=99; /$ This may not be what is desired. /$ This is the right way to do calculation using masks xx=(xx.eq.5)*99+(xx.ne.5)*xx; call print(x,xx); ENDIF End of an IF( )THEN structure endif; Must be the end of an if( )then; structure. Example: if(x.eq.0.0)then; y=10.; v=y*dsin(q); endif; PROGRAM Beginning of a program, program somename; begins a program. Programs use global variables. Example: program doit; call loaddata; call olsq(x,y :print); call graph(%res); return; end; RETURN( ) Returns the result of a function. return(result); Next to last statement in a user function. Example: function tt(y); t=sum(y); return(t); end; RETURN Next to last statement before end. return; must be the next to last statement in a PROGRAM, SUBROUTINE or FUNCTION. For a function the form return(value); is used. Example: program doit; call loaddata; call olsq(x,y :print); call graph(%res); return; end; SOLVE Solve a recursive system of equations. solve( ); solves a recursive system. solve(vv=x(t) :range 1 10); Solves recursively an expression. If formulas are involved, the form is solve(vv=x(t) :range 1 10 :block x); Note: :range i j is OK :range i, j is OK :range i, norows(gasout) is ok :range (i+6), j is not ok For more detail see extensive example below: b34sexec matrix; /$ Unlike RATS, SOLVE and FORMULA statements can /$ refer to themselves recursively n=1000; v=1.0; ar1=array(n:)+missing(); ar1(1)=99.+rn(v); solve(ar1=ar1(t-1)+rn(v):range 2 n); call graph(ar1); call tabulate(ar1); b34srun; SUBROUTINE Beginning of subroutine. subroutine somename(args); is the start of a subroutine. Subroutines use local variables. Example: subroutine test(x); call print('The mean of x is',mean(x)); return; end; WHERE( ) Starts a where structure. Note: The WHERE statement operates on non scalar objects while the IF statement operates on scalar objects. where(x.gt.0.0)y=x; Assuming x exists. where elements of x > 0.0 y=x, otherwise y=oldvalue. Variables x and y must be the same structure. The "mask" is done at the copy step. Commands of the form where(s.gt.0.0)y=dsqrt(s); will not work as intended since the dsqrt(s) is done BEFORE the mask is applied. Only simple replacement is allowed. Statements such as where(x.ge.0.0)z(,1)=q; are not allowed since z(,1) is a temp variable. The statements where(x.ge.0.0)x=missing(); where(x.ge.10.)x=dqsrt(20.); are allowed but care must be used. Since x and missing() and dsqrt(20.) are not the same structure, here if ( ) is "false" x is set = 0.0. For example, in the first statement where x ge 0.0, x is set to missing(). If x lt 0, x is set to 0.0. The second statement sets x to dsqrt(20.) if x was ge 10. Otherwise x is set to 0.0. Note that depending on whether x existed or not the statement does not return the same vector. This is an explicit design decision. Like Speakeasy, where( ) returns the existing value if the logical statement is false and zero if the variable did not exist! Assume x does not exist. Given where(x.ge.0.0)newx=missing(); Here if x ge 0.0, then newx = missing(), otherwise newx = 0.0. Logical masks are an alternative to some of the where capability. For example: x=a.gt.y; Here where a is > y, x = 1.0, x = 0.0 otherwise. Example: /; Here for the first where( ) the two objects /; across the equals sign are not the same structure /; If the ( ) is false x2bad resolves to 0.0 whether or /; not it existed prior to the where( ) being found. /; /; The second where( ) has objects the same structure across /; the =. Both objects exist. Here the old x value is /; maintained. The logic here is /; test = x*(x.ne.y)+dummy*(x.eq.y) b34sexec matrix; x=array(:1,-2,3,-4,5,-6,7,-8,9,-10); y=array(:0,-2,1,-4,6,-6,2,-8,5,-10); x2bad=x; x2good=x; dummy=array(norows(x):)+ -9999.; where(x.eq.y)x2bad =-9999.; where(x.eq.y)x2good =dummy; test = (x*(x.ne.y))+ (dummy*(x.eq.y)); call tabulate(x,y,dummy,x2bad,x2good,test); b34srun; Edited Results are: => X=ARRAY(:1,-2,3,-4,5,-6,7,-8,9,-10)$ => Y=ARRAY(:0,-2,1,-4,6,-6,2,-8,5,-10)$ => X2BAD=X$ => X2GOOD=X$ => DUMMY=ARRAY(NOROWS(X):)+ -9999.$ => WHERE(X.EQ.Y)X2BAD =-9999.$ => WHERE(X.EQ.Y)X2GOOD =DUMMY$ => TEST = (X*(X.NE.Y))+ (DUMMY*(X.EQ.Y))$ => CALL TABULATE(X,Y,DUMMY,X2BAD,X2GOOD,TEST)$ Obs X Y DUMMY X2BAD X2GOOD TEST 1 1.000 0.000 -9999. 0.000 1.000 1.000 2 -2.000 -2.000 -9999. -9999. -9999. -9999. 3 3.000 1.000 -9999. 0.000 3.000 3.000 4 -4.000 -4.000 -9999. -9999. -9999. -9999. 5 5.000 6.000 -9999. 0.000 5.000 5.000 6 -6.000 -6.000 -9999. -9999. -9999. -9999. 7 7.000 2.000 -9999. 0.000 7.000 7.000 8 -8.000 -8.000 -9999. -9999. -9999. -9999. 9 9.000 5.000 -9999. 0.000 9.000 9.000 10 -10.00 -10.00 -9999. -9999. -9999. -9999. *********************************************** Documentation of built-in commands called by CALL command. For futher examples, see problems in matrix.mac ACEFIT Alternating Conditional Expectation Model Estimation call acefit(y x[orderable] z[orderable]{2} :options); Implements the ace algorithm under matrix to provide estimation of ACE (alternating condition expectation) models following work by Breiman, L and Friedman J. (1985). Chapter 7 of Hastie and Tibshtiani (1990) provides a good reference. For another approach see the GAMFIT command. The ACEFIT command optionally can use the AVAS (Additivity and Variance Stabilization) approach. As noted in Hastie-Tibshirani (1990 page 193) "The AVAS procedure seeks transformations that achieve additivity and a homogenious variance and is more directed towards regresison problems than ACE." While AVAS may be more desirable for regresison models, its theoretical support is not as strong as for the ACE approach. - Faraway, Julian. "Extending the Linear Model with R" 2006, New York: Chapman & Hall/CRC - Hastie, T., R. Tibshirani and J. Friedman "The Elements of Statistical Learning: Data Mining, Inference, and Prediction." 2001 New York: Springer. - Breiman, L. and Friedman J. "Estimating Optimal Transformations for Multiple Regression and Correlation (with discussion),"Journal of American Statistical Association, 80, (1985) 580-619. - Hastie, T. J. and Tibshirani, R. J. (1990) "Generalized Additive Models," New York: Chapman and Hall Recent papers in the area include: - Fan, Jianqing and Jiancheng Jiang. "Nonparametric Inference for Additive Models," Journal of the American Statistical Association, 100, # 471 (2005) 890-907. where the generalized likelighood ratio test (GLR) is discussed. This test has not been implemented but can be implemented by interested users. Assume y = f(x1,...,x2) gamfit transforms the independent variables while acefit transforms both independent and dependent variables. Hastie and Tibshirani make the point that a model of the form y=exp(x1+x2**2)e cannot be estimated in additive form by gamfit but a simple additive model can be found that describes log(y). For example log(y)=x1+x2**2 + ln(e) Acefit procedure: Assume G(y) is the transformed dependent variable. The ACE procedure forces var{G(y)} = 1 and minimizes E{G(y)-f(X)}**2 subject to var{G(y)}=1. Assume one input variable x. yhat = ((G(y))**-1){a+sum(f(x))} where f(x]=) is the transformation for the input variable x and the transformation for the dependent variable is assumed to be invertable.. Logic of ACE: 1. For fixed G, the minimizing f is f(x)=E{G(y)|x} 2. For the minimizing f in step 1, the minimizing G is G(y)=E(f(x)|Y}. In this step new_(y)=estimated_g(y)/var{estimated_g(y)}**.5 Steps 1 and 2 alternate after assuming G(y) so that it has unit variance at each step to avoid a trivial zero solution. In words, the ace algorithm alternates between smoothing g(y) on x to get a new f(x), and f(x) on y to get a new g(y), until the mean-squared error does not change. Since the B34S implementation of ace saves all yhat vectors and residual vectors, the user can select the best model to use. The program call ace_ols; will save %best as a integer for the best model which allows the user to extract the transformed y and x matrix data and the %yhat vector for this model. Graphing of the transformed data can provide infcormation concerning nonlinearity in the model. The below listed examples show one criteria using the sum of squares of the residuals. However the user may want to weight the later performance more than the earlier performance of the model. The code to select best ace model (shown below) is placed after the call acefit command: /; Get best model call ace_ols; ace_res=%res(,%best); If the OLS model is not desired the call call ace_ols2(%k,%ns,%ssres,%ty,%tx,%best,%ybest,%xbest, iprint); can be used. To test the y variable transformation check to see if is_1= variancevariance(%besty)* (dfloat(norows(%besty)-1)/dfloat(norows(%besty))); is 1.0. Logic of AVAS In theory AVAS is more suitable for regession models in that it involves an asymptotic variance stabilizing transformation. For more detail see Stokes (200x) and Faraway (2006, 244-246). %tx and %ty can be studeied ands compared between the two procedures. Variables created %res - Residuals saved for last %nob %ns matrix %y - Y variable %yhat - Predicted y saved as a %nob,%ns array matrix %yvar - Y variable name %names - Names in Model %lag - Lag %vartype - Variable type %dist - Error Distribution %nob - Effective number of observations. %k - Number of right hand side variables. %ns - Number of passes. %rsq - Vector of R**2 for ns runs. This is the r**2 for the %ty series. The R**2 for y is given in the printout, if requested, but can be calculated from 1-(%ssres / variance(y)*dfloat(norows(y)-1)) %ssres - Vector of residual sum of squares. %best=imin(%ssres); %ty(n,ns) - Transformed y. %tx(n,p,ns) - Transformed x. [ ] Specifications. Allowed values are: order e => usual case. This is default for y and x. circular => periodic in range 0.0 to 1.0 with period 1. monotone => Transform must be monotone linear => transform is linear cat => not orderable. If the transformation is supplied, then the corresponding codes are fix_order fix_circular fix_monotone fix_linear fix_cat Examples: Call acefit(y x1 x2[order] x3[order]{1} x4[order]{1 to 6} :print); Note: If lags are present then based on minimum lag the following is saved. This will not occur if future right hand side variables are present and :holdout is not present. %xfobs - Observation number. %xfuture - Same as %x but for out of sample data that is available. Note: while x1 is allowed x1{1} is not since [ ] is missing. Options supported :print Show output :itprint Show Iterations :avas Use AVAS approach to get transformations. This model does not forecast at this time. It is more a diagnostic tool. :span r1 Sets span. Default = 0.0. Set as a fraction of observations in window. 0.0 => automatic (variable) span selection. For small samples (n < 40) or if there are substantial serial correlations between obserations close in x - value, then a prespecified fixed span smoother (span > 0) should be used. Reasonable span values are 0.3 to 0.5. :alpha r2 Sets alpha. Default = 0.0. Controls high frequency (small span) penality used with automatic span selection (base tone control). alpha.le.0.0 or alpha.gt.10.0 => no effect.) :ns ns Number of solutions. Default = 3. :ty ty Supplied y transformation. ty(n,ns) :tx tx(n,p,ns) Supplied x transform :maxit i1 Maximum number of iterations. Default=20 :nterm i2 Maximum number of terminal iterations. Default = 3 :tol delrsq Sets termination threshold. Default=.1e-3. Iteration stops when R**2 changes less than delrsq in i2 iterations. :holdout n Sets number of observations to hold out :sample mask Specifies a mask real*8 variable that if = 0.0 drops that observation. Note: :sample cannot be used with :holdout. Unless the mask is the number of obs after any lags, an error message will be generated. The sample variable must be used with great caution when there are lags. A much better choice is the :holdout option. :weight w Set weight for each obs. Series w must have same # of observations as left hand variable. :spans r3 Sets a three element array. Spans values are for the three running linear smoothers. spans(1) : tweeter span. spans(2) : midrange span. spans(3) : woofer span. Default => array(:.05,.2,.5) This parameter should not be adjusted under normal circumstances. :savex Saves the X matrix in %X. call acefit(y x[order]{1 to 6} z[cat] :print); The model specification involves specification of the type of variable and optionally a lag or lags. The model specification allows the lags to be set in the command. Only vectors can be supplied in this release. If no[ ] is supplied, [order] is assumed. The specification call acefit(y y[order]{1} x[order]{0 to 3} z[order]{1} )$ is the same as call acefit(y y[order]{1} x[order] x[order]{1} x[order]{2} x[order]{3} z[order]{1})$ :getmodel 'filename' Saves model discription variables If getmodel is found, no estimation is performed. :savemodel 'filename' The default name is 'acemod.psv' :forecast xmatrix => Allows users to supply observations of the right hand side variables outside the sample period so that forecasts can be calculated. The same number of observations must be supplied for all right hand series. Due to the way that splines are calculated, it is imperative that any values of the x variables NOT lie outside the ranges of the original data. The user must have save the model workspace if :forecast is the only option. The forecast sentence produces the %fore variable and the %foreobs variable. As of this time :forecast works only for ACE models. :fmodel integer Sets model number to forecast. Default=1 Variables saved if :savemodel is set: %X, %K, %YVAR, %NAMES, %LAG, %NS, %TY, %TX, %Y, %NOB, %RSQ, %RES, %YHAT, %LF, %M, %Z, %F, %T Many of these values have use only in forecasting. Important supplementary routines in staging routine program ace_ols; /; /; Extracts "best" ACE Model - saves transformed y in /; %ybest and saves transformed x in %xbest %best = index /; /; Obtains %k, %ssres %ty and %tx from last acefit model /; /; Experimental version 19 February 2006 /; /; For reference and motivation see Faraway (2006) P242 /; Can be called as a program or as a subroutine using iprint=1; call aceols2(%k,%ns,%ssres,%yhat,%ty,%tx,%best,%ybest, %xbest, iprint,%ace,%y); call print(' ':); call print('Test sum of squares ',sumsq(%y-%ybest):); return; end; subroutine aceols2(k,ns,ssres,yhat,ty,tx,best,ybest,xbest, iprint,ace,y); /; /; Extracts best ACE model data and optionally shows /; implied OLS models /; /; Let y* = theta(y) /; /; 1. y* = f( smoothed x) /; 2. y* = f( smoothed x without constant) /; 3. y = f( smoothed x ) /; /; Y can be forecasted two ways: /; -using inverse of y filter /; -using equation 3. /; /; Can be called by ACE_OLS or directly /; /; k => # right hand side variables /; ns => # of ACE models attempted /; ssres => array with ns elements showing residual sum /; of squares /; of ACE model /; yhat => nob,ns 2d object containing y hat /; ty => transformed y /; tx => nob,k,ns of nob * (k*ns) object containing /; transformed x /; best => # of ACE model that has minimum e'e /; ybest => nob element transformed y of "best" model /; xbest => nob,k matrix of transfrmed x elements of /; "best" model /; iprint => =0 only calculate, NE 0 => print OLS models /; of ybest = f(xbest) /; ace => Set by call acefit. use %ACE /; /; Experimental version 19 February 2006 /; /; For reference and motivation see Faraway (2006) P242 /; program ace_plot; /; /; Plot ACE Transformations for BEST Model /; subroutine aceplot2(y,ybest,xbest,xorig,names,lag,ace); /; /; Show Ace Transforms in Graphs /; See Faraway (2006) page 243 /; /; y => Usually %y /; ybest => Best ytransform /; xbest => best xtransform /; xorig => %x from :savex /; names => %names /; lag => %lag /; ace => %ACE Examples of using supplementary routines. Shows GAM vs ACE model. /; /; ACE Plots /; b34sexec options ginclude('b34sdata.mac') member(b_ozone); b34srun; b34sexec matrix; call loaddata; call echooff; call load(ace_ols :staging); call load(ace_plot :staging); iholdout=10; call gamfit(ozone vh[predictor,3] wind[predictor,3] humidity[predictor,3] temp[predictor,3] ibh[predictor,3] dbg[predictor,3] ibt[predictor,3] vis[predictor,3] :holdout iholdout :print); call acefit(ozone vh[order ] wind[order ] humidity[order ] temp[order ] ibh[order ] dbg[order ] ibt[order ] vis[order ] :holdout iholdout :ns 4 :savex :print); call ace_ols; call ace_plot; b34srun; Examples: b34sexec options ginclude('b34sdata.mac') member(gam); b34srun; b34sexec options noheader; b34srun; b34sexec matrix; call loaddata; call echooff; call acefit(y[cat] age[order] start_v[order ] numvert[order] :print); call gamfit(y age[predictor,3] start_v[predictor,3] numvert[predictor,3] :link logit :dist gauss :maxit index(2000,1500) :tol array(:.1d-13,.1d-13)); b34srun; /; /; ACEFIT Problem showing plots for various NS values /; b34sexec options ginclude('b34sdata.mac') member(gam_3); b34srun; b34sexec options noheader; b34srun; b34sexec matrix; call loaddata; call olsq(cpeptide age bdeficit : print); call echooff; call acefit( cpeptide[order ] age[order] bdeficit[order] :maxit 20 :nterm 10 :ns 2 :tol .1e-8 :print); call names(all); call tabulate(%rsq,%ssres); call print(%y); call print(%yhat,%res); do i=1,%ns; call graph(%res(,i)); enddo; b34srun; /; /; Experimental AVAS Option /; b34sexec options ginclude('b34sdata.mac') member(gam_3); b34srun; b34sexec options noheader; b34srun; b34sexec matrix; call loaddata; call olsq(cpeptide age bdeficit : print); call echooff; call acefit( cpeptide[order ] age[order] bdeficit[order] :avas :maxit 20 :nterm 10 :tol .1e-8 :print); call names(all); call tabulate(%rsq,%ssres); call print(%y); call print(%yhat,%res); do i=1,%ns; call graph(%res(,i)); enddo; b34srun; /; /; ACEFIT Vs GAMFIT on Gas Data /; b34sexec options ginclude('b34sdata.mac') member(gas); b34srun; b34sexec options noheader; b34srun; b34sexec matrix; call loaddata; call load(ace_ols :staging); maxlag=6; call olsq(gasout gasout{1 to maxlag} gasin{1 to maxlag} :print); ols_res=%res; call acefit(gasout gasout[order]{1 to maxlag} gasin[order]{1 to maxlag} :print); call ace_ols; ace_res=%res(,%best); call gamfit(gasout gasout[predictor,3]{1 to maxlag} gasin[predictor,4]{1 to maxlag} :print); gam_res=%res; call marspline(gasout gasin{1 to 4} gasout{1 to maxlag} :nk 40 :mi 2 :print); mars_res=%res; call graph(ols_res,gam_res,ace_res,mars_res :nolabel); Models=c8array(:'OLS','ACE','GAM','MARS'); fit =array(:sumsq(ols_res),sumsq(ace_res), sumsq(gam_res),sumsq(mars_res)); call tabulate(models,fit); b34srun; Forecasting ACE vs MARS /; /; ACE validation of Forecasting - /; Feeds data back in to prove forecast /; Out of sample Forecast done /; /; Data set used was studied in Faraway (2006) page 241 b34sexec options ginclude('b34sdata.mac') member(b_ozone); b34srun; b34sexec matrix; call loaddata; call echooff; iholdout=80; idoace =1; iacefore=1; idomars=1; /; ACE +++++++++++++++++++++++++++++++++++++++++++++++++++ if(idoace.ne.0)then; call acefit(ozone vh[order ] wind[order ] humidity[order ] temp[order ] ibh[order ] dbg[order ] ibt[order ] vis[order ] doy[order ] :savemodel :holdout iholdout :ns 1 :savex :print); %x_ace=%x; %yhatt=%yhat; /; First forecast the complete sample Check sum of squares if(iacefore.ne.0)then; call acefit(ozone vh[order ] wind[order ] humidity[order ] temp[order ] ibh[order ] dbg[order ] ibt[order ] vis[order ] doy[order ] :holdout iholdout :getmodel :forecast %x_ace :print); call print('%fore and %yhatt have to be the same':); call tabulate(%foreobs %fore,%yhatt,ozone); jj=integers(norows(ozone)-iholdout); %actual=ozone(jj); call print(' '); call print('Test Error Sum of squares both ways ':); call print('Error sum of squares using yhat ',