参考：http://140.123.5.6/deptfin/course_data/Gauss.htm

cls;

a = { 2.2 1.7 211221};   /* Generating 1*3 MATRIX */

b = {-3.2 3.9  43433434};

hc = a ~ b;        /* 横向连接矩阵 */

print hc;   /* 计算机屏幕 显示数据 */

vc = a|b;   /* 数据纵向 连接 */

print vc;

——————————————————————————————————————————————

new；

cls;

proc (1) = demean(x);   /*  这个相当于一个SAS宏语句，或者MATLAB,R 的函数*/

local nobs;

local xbar;

xbar = meanc(x);

retp(x - xbar);

endp;

load path=d: ;      /*  设置工作路径 ，然后读入 TXT 数据（纯数据）*/

load x[1200,6]=aaa.txt;

print x;

data = x[2:10,.];

print data;

consdm = demean(data[.,2]); /*  运行前面的 宏代码 */

print consdm; data[.,2];

e.g. ：

proc(1)=detpd(X) ; /* (1) is the number of returned outputs */

local Y ; /* note that you have to write down the local var you are using in your proc

this can be done in several steps as here or using: local Y, X ; */

local Z ;

Y=X'*X ; /* may be faster to use =MOMENT(X,0) */

Z=det(Y) ;

retp(Z) ; /* specifies the output of the proc */

endp ; /* note the syntax proc(#arg)=func() ; endp ;*/

X=ones(5,2) ; /* let’s test the process */

d0=detpd(X) ;

print "the det of X'X is" ;; print d0 ;

————————————  储存 与 载入 GAUSS 格式 数据 ————————————————

1. ASCII 数据的 调入 GAUSS 系统：

方式一： load x[N, K]=datafile.asc;

方式二：load x[]=datafile.asc;

       x=reshape(x,N,K);

2.将调入 GAUSS 系统的数据保存为 GAUSS 格式：

cls;

Let x={0 2,2.0000001 4.00001} ;

Let xnames={"X1","X2"} ;

Export=saved(x,"D:\example",xnames) ;  /* 这样就保存为GAUSS格式了 example ，外形还是.dat*/

if export==0 ;

"the export is unsuccessful" ;

else ;

"the export is successful" ;

endif ;

3. 载入之前 或者 他人的GAUSS 格式数据：

方式一：

new ;

CLOSEALL ;  /* this closes all datasets previously opened */

/* the double backward slashes are needed to specify the path within the quotes */

/* creates a string macro filename */

/* the name of the storage file is lqdata */

/* gauss stores the data into 2 parts .dat=matrix, .dht=var names */

open input=d:\example ; /* opens the dataset */

data2=READR(input,ROWSF(input)) ; /* reads it into some matrix called data */

/* rowsf(input) specify we want the S observations. */

/* It can be changed for 10, 100 etc. */

vnames=getname("D:\example") ; /* takes the names of the variables */

print data2[.,1] ;  /*显示 调入 工作空间 的数据*/

print $vnames;     /*显示调入工作空间的数据——变量 名称*/

方式二：

yy=loadd("D:\example");

yy;

/* GAUSS is not case-sensitive. */

_________________________   矩阵的运算  ——————————————————

a) Fixed values matrices

A=zeros(10,3) ; /* creates a matrix of 0s 10 rows * 3 columns */

I=eye(10) ; /* creates an identity matrix 10*10 */

B=ones(5,12) ; /* creates a matrix of 1s 5 rows * 12 columns */

C=seqa(start,inc,n) ; /* creates an additive column vector sequence of increment

inc, starting at C1=start, C2=start+inc… and containing n values */

D=seqm(start,inc,n) ;/* creates a multiplicative vector sequence as seqa */

The user can also create a matrix using constants with the LET command:

LET x=1 2 3 4 ; /* column vector */

LET x=1,2,3,4 ; /* as before */

LET x=1 2,3,4 ; /* as before! */

LET x={1 2 3 4} ; /* row vector */

LET M={1 2, 3 4} ; /* 2*2 matrix with 1st row 1 2 */

LET M[2,3]=1 2 3 4 5 6 ; /* 2*3 matrix with 1st row 1 2 3*/

LET M[2,3]=1, 2, 3, 4, 5, 6 ; /* 2*3 matrix with 1st row 1 2 3*/

LET M[2,3]= 6 ; /* 2*3 matrix with all entries=6 */

b) “Random” matrices

N=RNDN(10,3) ; /* creates a 10*3 matrix of independent pseudo random draws of a

N(0,1) */

S=1 ; N1=RNDNS(10,3,S) ; /* creates a 10*3 matrix of independent pseudo random

draws of a N(0,1) with a given seed S=1*/

U=UNDN(10,3) ; /* creates a 10*3 matrix of independent pseudo random draws of a

Uniform (0,1) */

S=12 ; U1=UNDNS(10,3,S) ; /* creates a 10*3 matrix of independent pseudo random

draws of a Uniform (0,1) with a given seed S=12*/

LET V={1 0.99,0.99 1} ;

Sl=chol(V)' ; /* Sl is Cholesky lower triangular V1/2 in some textbooks */

print V ;

print Sl ;

Seed=101 ;

U=rndns(2,100,Seed) ;

B = Sl*U ;

G=U|B ; /* vertical concatenation for comparison */

print /* U|B */ ; /* the second series of draws B1, B2 has values of B1 very close to

B2 at each draw, while U1 and U2 are unrelated */

print G ;

3.2. Operations on matrices5

a) Referencing matrices

This works in the same way for rows and columns:

g=mat[r1:r2,c1:c2] ; /* gives you rows r1 to r2 and col c1 to c2 of the matrix mat */

g=mat[r,c1:c2] ; /* gives you row r and columns c1 to c2 of the matrix mat */

g=mat[r,c1 c2 c3] ; /* gives you row r and columns c1, c2, c3 of the matrix mat */

g=mat[.,c1:c2] ; /* gives you all the rows and columns r1 to r2 of the matrix mat */

b) Simple operations

C=A*B ; /* product */

C=A’ ; /* transpose */

Rk: Computing Z=X’X seems quite useful but GAUSS has an in built routines that is

much faster:

Z=MOMENT(X,0) ; /* returns X’X ignoring missing values (1 or 2 instead of 0

delete them) */

C=inv(A) ; /* inverse */

C=invpd(A) ; /* inverse of a pos def matrix eg. X’X */

D=diag(A) ; /* diagonal of A */

D=det(A) ; /* det of A */

D=rank(A) ; /* rank of A */

C=cols(A) ; /* cols of A */

R=rows(A) ; /* rows of A */

E=eig(A) ; /* eigenvalues of A */

Pi is pi approx. 3.1415927.

^ is exponentiation

! is factorial

% is modulo division

C=A|B ; /* vertical concatenation of A and B */

C=A~B ; /* horizontal concatenation of A and B */

C=meanc(A) ; /* returns a column vector of the column mean of A */

Rk: Maxc(), Minc(),Sumc(), stdc() and median() do the same for the related

operations.

C=A.*.B ; /* Kronecker product of A and B */

C=Vec(A) ; /* returns the column vector of the stacked columns below one another */

C=Vech(A) ; /* takes the lower triangular part of a symmetric matrix in a column

vector */

Su=Chol(V) ; /* returns the upper triangular part of the cholesky decomp V=Su’Su */

It is also possible to sort a matrix with respect to a particular column (sortc for

numeric values or sortcc for character values) or a group of columns (sortmc). The

following illustrates the ranking of the values in GAUSS:

Let X={2 2 5 "E", . 1 6 "b", 2 1 3 "c", 3 2 7 "B", 4 3 8 "D"} ;

let mask=1 1 1 0 ; /* we define the columns of the matrix as num (1) or char (0) */

printfmt(X, mask') ; /* we print the matrix of char and num values */

X1=sortc(X,1) ; /* we sort wrt 1st column which contains numeric values */

X2=sortcc(X,4) ; /* we sort wrt 4th column which contains char. values */

printfmt(X1, mask') ; /* missing values are considered as arbitrarily small */

printfmt(X2, mask') ; /* GAUSS considers All capitalized letters as smaller as the

un capitalized ones */

let sorter=1 2 ; /* we plan to sort the matrix wrt columns 1 2 of num. values */

X3=sortmc(X,sorter) ; printfmt(X3, mask') ;

let sorter=1 -4 ; /* we plan to sort the matrix wrt columns 1 4 of num. values and

char. Values so we had a – before the 4 */

X4=sortmc(X,sorter) ; printfmt(X4, mask') ;

c) Operations element by element

C=sqrt(A) ; /* returns sqrt(aij) other functions are ln, log, exp, cos, sin */

C=A^2 ; /* returns aij^2 */

C=A./B ; /* returns aij/bij it works also with B a vector and gives a matrix aij/bi */

C=A.*B ; /* returns aij*bij it works also with B a vector and gives a matrix aij*bi */

Rk In general you can use the usual symbol but with a dot to do element by element

operation. Eg. A.^B gives aij^bij.

C=A-3 ; /* returns aij-3 */

C=A+3 ; /* returns aij+3 */

C=A*3 ; /* returns aij*3 */

C=A/3 ; /* returns aij/3 */

d) Relational and logical operators

== or EQ for equality ( this is different from A=5 where A is assigned value 5)

/= or NEQ for not equal

> or GT for greater than

>= or GE for greater or equal than

< or LT for lower than

<= or LE for lower or equal than

GAUSS returns 1 if a statement is TRUE and 0 if FALSE. You can combine

statement using logical operators:

NOT Statement 1

Statement1 AND statement2

Statement1 OR statement2

Statement1 XOR statement2 (one, or the other, but not both)

Eg.

LET V={1 0.8,0.8 1} ;

LET V1={1 0.7,0.8 1} ;

C=V.==V1 ; /* this check element by element the == of the 2 matrices */

Print C ;

This returns:

1.0000000 0.0000000

1.0000000 1.0000000

In numerical computations, it is useful to be more tolerant on the meaning of ==.

There is always some elements of rounding. You can set yourself the comparison with

some tolerance level or use the fuzzy comparison operators of GAUSS. This

corresponds to the usual operators with a prefix F (eg EQ is now FEQ). This can also

be use for element by element comparison using the prefix DOTF (eg. DOTFEQ).

Eg.

A=1*10^(-15) ; B=1*10^(-16) ; TOL=1*10^(-15) ;

C=A==B ;

CT=A GE B-TOL AND A LE B+TOL ; /* this uses the logical operator AND !! */

_fcmptol=10^(-15) ; /* set the tolerance limit of gauss to 10^-15 this is the default

value for the fuzzy comparison operators */

CTG=FEQ(A,B) ;

Print C ;; Print CT ;; Print CTG ;

Returns :

0.00000000 1.0000000 1.0000000

>> A={1 2 3,1 4 5};

>> A;

      1.0000000        2.0000000        3.0000000

      1.0000000        4.0000000        5.0000000

>> B={1 4 5, 4 8 0};

>> B;

      1.0000000        4.0000000        5.0000000

      4.0000000        8.0000000       0.00000000

>> A*B';

      24.000000        20.000000

      42.000000        36.000000

>>A.*.B;

    1.0      4.0      5.0      2.0      8.0      10.      3.0      12.      15.

    4.0      8.0     0.00      8.0      16.     0.00      12.      24.     0.00

    1.0      4.0      5.0      4.0      16.      20.      5.0      20.      25.

    4.0      8.0     0.00      16.      32.     0.00      20.      40.     0.00

>> A.*B;

      1.0000000        8.0000000        15.000000

      4.0000000        32.000000       0.00000000

>> c=seqa(2,0.25,10) ;  /* 连续 相加 数列*/

>> c;

      2.0000000

      2.2500000

      2.5000000

      2.7500000

      3.0000000

      3.2500000

      3.5000000

      3.7500000

      4.0000000

      4.2500000

>> D=seqm(2,0.25,10) ;   /* 连续 相乘 数列*/

>> d;

      2.0000000

     0.50000000

     0.12500000

    0.031250000

   0.0078125000

   0.0019531250

  0.00048828125

  0.00012207031

 3.0517578e-005

 7.6293945e-006

———   Procedures ：4.1. Writing and using a procedure (PROC ENDP) ————

A procedure is a sub-program which will help you doing repeated tasks inside the

same program. It has a number of inputs or arguments and can return several outputs.

A simple procedure may have only one argument. For example, we can create a

procedure that multiplies X’ by X and return the determinant of this matrix.

proc(1)=detpd(X) ; /* (1) is the number of returned outputs */

local Y ; /* note that you have to write down the local var you are using in your proc

this can be done in several steps as here or using: local Y, X ; */

local Z ;

Y=X'*X ; /* may be faster to use =MOMENT(X,0) */

Z=det(Y) ;

retp(Z) ; /* specifies the output of the proc */

endp ; /* note the syntax proc(#arg)=func() ; endp ;*/

X=ones(5,2) ; /* let’s test the process */

d0=detpd(X) ;

print "the det of X'X is" ;; print d0 ;

———————— 使用 GAUSS自带的 包 Some existing procedures (OLS, DSTAT) —————

new;

cls;

x=RNDN(1000000,1);

y=RNDN(1000000,1);

{vnam,m,b,stb,var,stddev,sig,cx,rsq,res,dwstat}=ols(0,y,x) ; /* see help OLS */

u=rndn(1000,4);

{vnam,mean,var,std,min,max,valid,mis}=dstat(0,u);

xy=x~y;

xy1=xy[1:1000,1];

{vnam,mean,var,std,min,max,valid,mis}=dstat(0,xy1);

________________________   Bootstrap 标准误回归 OLS   _____________________________

/*

**Bootstrap OLS

*/

@ Data generation under Weak Ideal Conditions @

/*

** Regressors are stochastic.

** The errors are different across different data sets.

** Errors are non-normal: Chi-square(2)-2

*/

new;

cls;

@ open output file @

   output file = bootols_cross.out reset ;

@  Model: y(t) = beta(1) + beta(2)*x(t) + e @

seed   = 1;

beta2  = 0.5;

beta1 = 0.5;

tt      = 100; @ # of observations @

kk      = 2;  @ # of betas @

nboot   = 1000;

@ Generating x and y  @

x = rndns(tt,1,seed)^1;

y = beta1 + beta2*x + (rndns(tt,1,seed)^2+rndns(tt,1,seed)^2-2);

@OLS using yy and xx  @

yy = y;

xx = ones(tt,1)~x;

olsb  = invpd(xx'xx)*(xx'yy);

olse  = yy - xx*olsb ;

s2 = (olse'olse)/(tt-kk);

v  = s2*invpd(xx'xx);

se = sqrt(diag(v));

@ testing H_o: beta2 = true beta 2 @

asympw = ( (olsb[2]-beta2)/se[2] )^2 ;

asympp = cdfchic(asympw,1);

asympw1 = ( (olsb[2]-1)/se[2] )^2 ;

asympp1 = cdfchic(asympw1,1);

"Asymptotic Wald test result for H_o: beta2 = true beta2 (" beta2 ")";

"  Wald Stat,  p-val:" asympw asympp ;

"";

"Asymptotic Wald test result for H_o: beta2 = 1";

"  Wald Stat,  p-val:" asympw1 asympp1 ;

"";

@ Boot strap tests @

bootw = zeros(nboot,1) ;

z = y~x;

i = 1 ;

do while i <= nboot ;

   cboot = ceil( tt*rndus(tt,1,seed) ) ; @ choosing bootstrap data points @

   bz    = z[cboot,.] ;

   byy = bz[.,1];

   bxx = ones(tt,1)~bz[.,2];

   @Bootstrap OLS using yy and xx  @

   bootb  = invpd(bxx'bxx)*(bxx'byy);

   boote  = byy - bxx*bootb ;

   bs2 = (boote'boote)/(tt-kk);

   bv  = bs2*invpd(bxx'bxx);

   bse = sqrt(diag(bv));

   bootw[i] = ( (bootb[2]-olsb[2])/bse[2] )^2 ;

i = i + 1 ;

endo ;

@ BOOTSTRAP CRITICAL VALUES @

   bootwt = sortc(bootw,1) ;

   tc     = nboot           ;

   bootc  = bootwt[ceil(tc*0.99)]|bootwt[ceil(tc*0.95)]|bootwt[ceil(tc*0.9)];

"" ;

"BOOTSTRAP WALD TEST FOR H_o: beta2 = true (" beta2 ")" ;

"";

"       STATISTICS       C-VAL-1%         C-VAL-5%         C-VAL-10% ";

asympw bootc' ;

"" ;

"BOOTSTRAP WALD TEST FOR H_o: beta2 = 1";

"";

"       STATISTICS       C-VAL-1%         C-VAL-5%         C-VAL-10% ";

asympw1 bootc' ;

output off;