BookmarkSubscribeRSS Feed
tammytt
Fluorite | Level 6

Hi, 

I am running a bunch of code which is published by a researcher on the sas of wrds cloud and get the error of the following. 

I am totally new to sas and have no idea where the problem could come from. Any hints will be really appreciated!

 

210        *____________________________________________________________________________________________________
211        
212        	This proc loops over all samples, calling the estimation routines.
213        _____________________________________________________________________________________________________;
214        options nonotes;
215        proc iml;
216        	
216      !  start main;
217        	
217      !  call streaminit(1234);
218        	
218      !  reset storage=this.imlstor;
218      !                             	*	This contains the subroutines built in RollGibbsLibrary02.sas;
219        	
219      !  load;
220        
221        	
221      !  reset printadv=1 log;
222        	
222      !  use temp.pys;
223        	
223      !  read all var {permno year kSample} where(nTradeDays>=60) into sample [colname=colSample];
223      !                                                                                             *my correction;
224        	
224      !  nSamples = nrow(sample);
225        	
225      !  print 'nSamples=' nSamples;
226        	
226      !  permno = sample[,1];
227        	
227      !  year = sample[,2];
228        	
228      !  kSample = sample[,3];
229        	
229      !  outSet = j(1,7,.);
230        	
230      !  varnames ={'permno','year','kSample','c','beta','varu','sdu'};
231        	
231      !  create this.gibbsOut from outSet [colname=varnames];
9 The SAS System                                                                               Sunday, August  8, 2021 02:08:00 AM

232        	
232      !  do iSample=1 to nSamples;
233        		
233      !   thisPermno = permno[iSample];
234        		
234      !   thisYear = year[iSample];
235        		
235      !   thisKSample = kSample[iSample];
236        		
236      !   if mod(iSample,500)=1 then do;
237        			
237      !    t = time();
238        			
238      !    ctime = putn(t,'time.');
239        			
239      !    print ctime iSample '/' nSamples ': ' thisPermno thisYear thisKSample;
240        		
240      !   end;
241        		
241      !   if thisPermno>=&startPermno & thisYear>=1990 then do;
242        			
242      !    use temp.dsf where (permno=thisPermno & year=thisYear & kSample=thisKSample);
243        			
243      !    read all var {p pm q} into x [colname=colx];
244        			
245        			
245      !    nSweeps = 1000;
246        			
246      !    regDraw = 1;
247        			
247      !    varuDraw = 1;
248        			
248      !    qDraw = 1;
249        			
249      !    nDrop = 200;
250        		
251        			
251      !    call RollGibbsBeta(parmOut, x[,1],x[,2],x[,3], nSweeps, regDraw, varuDraw, qDraw, 0,0,0,0);
252        			
253        			
253      !    p2 = parmOut[(nDrop+1):nSweeps,];
254        			
254      !    p2 = p2 || sqrt(p2[,3]);
255        			
255      !    pm = p2[+,]/(nSweeps-nDrop);
256        			
256      !    outset = thisPermno || thisYear || thisKSample || pm;
257        			*print outset;
258        			
258      !    setout this.gibbsOut;
259        			
259      !    append from outSet;
260        		
260      !   end;
261        	
261      !  end;
262        	
262      !  finish main;
263        run;
10 The SAS System                                                                              Sunday, August  8, 2021 02:08:00 AM

           nSamples

nSamples=    356031

ctime      iSample    nSamples    thisPermno  thisYear thisKSample

 2:19:25         1 /    356031 :       10001      1986           1

ctime      iSample    nSamples    thisPermno  thisYear thisKSample

 2:21:07       501 /    356031 :       10034      1986           1

ctime      iSample    nSamples    thisPermno  thisYear thisKSample

 2:22:58      1001 /    356031 :       10078      1987           1

ctime      iSample    nSamples    thisPermno  thisYear thisKSample

 2:25:16      1501 /    356031 :       10121      1988           1

ctime      iSample    nSamples    thisPermno  thisYear thisKSample

 2:27:09      2001 /    356031 :       10153      1966           1

ctime      iSample    nSamples    thisPermno  thisYear thisKSample

 2:28:56      2501 /    356031 :       10196      1950           1

ctime      iSample    nSamples    thisPermno  thisYear thisKSample

 2:30:55      3001 /    356031 :       10233      1938           1

ctime      iSample    nSamples    thisPermno  thisYear thisKSample

 2:33:04      3501 /    356031 :       10268      1950           1

ctime      iSample    nSamples    thisPermno  thisYear thisKSample

 2:35:17      4001 /    356031 :       10308      2014           1

ctime      iSample    nSamples    thisPermno  thisYear thisKSample

 2:37:34      4501 /    356031 :       10362      2002           1

ctime      iSample    nSamples    thisPermno  thisYear thisKSample

 2:39:37      5001 /    356031 :       10401      1930           1

ctime      iSample    nSamples    thisPermno  thisYear thisKSample

 2:41:34      5501 /    356031 :       10443      2017           1

ctime      iSample    nSamples    thisPermno  thisYear thisKSample

 2:42:51      6001 /    356031 :       10487      1967           1

ctime      iSample    nSamples    thisPermno  thisYear thisKSample

 2:44:33      6501 /    356031 :       10527      1995           1

ctime      iSample    nSamples    thisPermno  thisYear thisKSample

 2:46:41      7001 /    356031 :       10571      2002           1

ctime      iSample    nSamples    thisPermno  thisYear thisKSample

 2:48:17      7501 /    356031 :       10620      1940           1
11 The SAS System                                                                              Sunday, August  8, 2021 02:08:00 AM

ctime      iSample    nSamples    thisPermno  thisYear thisKSample

 2:50:41      8001 /    356031 :       10661      1998           1

ctime      iSample    nSamples    thisPermno  thisYear thisKSample

 2:52:46      8501 /    356031 :       10715      1988           1

ctime      iSample    nSamples    thisPermno  thisYear thisKSample

 2:53:56      9001 /    356031 :       10755      1994           1

ctime      iSample    nSamples    thisPermno  thisYear thisKSample

 2:56:07      9501 /    356031 :       10800      1996           1

ctime      iSample    nSamples    thisPermno  thisYear thisKSample

 2:58:27     10001 /    356031 :       10853      2012           1

ctime      iSample    nSamples    thisPermno  thisYear thisKSample

 3:00:25     10501 /    356031 :       10890      1956           1

ctime      iSample    nSamples    thisPermno  thisYear thisKSample

 3:02:51     11001 /    356031 :       10933      1995           1

ctime      iSample    nSamples    thisPermno  thisYear thisKSample

 3:04:47     11501 /    356031 :       10984      1993           1

ctime      iSample    nSamples    thisPermno  thisYear thisKSample

 3:06:39     12001 /    356031 :       11033      1963           1

ctime      iSample    nSamples    thisPermno  thisYear thisKSample

 3:08:40     12501 /    356031 :       11078      1994           1

ctime      iSample    nSamples    thisPermno  thisYear thisKSample

 3:10:12     13001 /    356031 :       11135      1989           1

ctime      iSample    nSamples    thisPermno  thisYear thisKSample

 3:12:07     13501 /    356031 :       11174      2014           1

ctime      iSample    nSamples    thisPermno  thisYear thisKSample

 3:14:20     14001 /    356031 :       11244      1951           1

ctime      iSample    nSamples    thisPermno  thisYear thisKSample

 3:16:16     14501 /    356031 :       11293      1988           1

ctime      iSample    nSamples    thisPermno  thisYear thisKSample

 3:17:46     15001 /    356031 :       11335      1997           2

ctime      iSample    nSamples    thisPermno  thisYear thisKSample

 3:19:40     15501 /    356031 :       11368      2018           1

ctime      iSample    nSamples    thisPermno  thisYear thisKSample

 3:22:04     16001 /    356031 :       11403      1994           1
12 The SAS System                                                                              Sunday, August  8, 2021 02:08:00 AM

ctime      iSample    nSamples    thisPermno  thisYear thisKSample

 3:23:58     16501 /    356031 :       11447      1975           1

ctime      iSample    nSamples    thisPermno  thisYear thisKSample

 3:25:49     17001 /    356031 :       11490      1995           1

ctime      iSample    nSamples    thisPermno  thisYear thisKSample

 3:27:43     17501 /    356031 :       11535      1936           1

ctime      iSample    nSamples    thisPermno  thisYear thisKSample

 3:29:55     18001 /    356031 :       11599      1993           3

ctime      iSample    nSamples    thisPermno  thisYear thisKSample

 3:32:18     18501 /    356031 :       11643      1991           1

ctime      iSample    nSamples    thisPermno  thisYear thisKSample

 3:34:23     19001 /    356031 :       11678      2001           1

ctime      iSample    nSamples    thisPermno  thisYear thisKSample

 3:36:34     19501 /    356031 :       11729      1996           1

ctime      iSample    nSamples    thisPermno  thisYear thisKSample

 3:38:31     20001 /    356031 :       11764      1988           1

ctime      iSample    nSamples    thisPermno  thisYear thisKSample

 3:41:03     20501 /    356031 :       11826      1930           1

ctime      iSample    nSamples    thisPermno  thisYear thisKSample

 3:43:01     21001 /    356031 :       11868      1995           1

ctime      iSample    nSamples    thisPermno  thisYear thisKSample

 3:45:34     21501 /    356031 :       11914      1951           1

ctime      iSample    nSamples    thisPermno  thisYear thisKSample

 3:48:02     22001 /    356031 :       11981      1926           1

ctime      iSample    nSamples    thisPermno  thisYear thisKSample

 3:50:21     22501 /    356031 :       12018      2012           1

ctime      iSample    nSamples    thisPermno  thisYear thisKSample

 3:52:11     23001 /    356031 :       12052      1994           2

ctime      iSample    nSamples    thisPermno  thisYear thisKSample

 3:54:45     23501 /    356031 :       12079      1931           1

ctime      iSample    nSamples    thisPermno  thisYear thisKSample

 3:56:25     24001 /    356031 :       12110      2002           1
ERROR: (execution) Module not loaded, operation not available.

 operation : RAND at offset  10 column  16
 operands  : *LIT1142

13 The SAS System                                                                              Sunday, August  8, 2021 02:08:00 AM

*LIT1142      1 row       1 col     (character, size 7)

 uniform

 statement : ASSIGN at offset  10 column   1
 traceback : module RANDSTDNORMT at offset  10 column   1
             module MVNRNDT at offset  12 column   2
             module ROLLGIBBSBETA at offset  39 column   4
             module MAIN at line 251 column 4

264        quit;
6 REPLIES 6
JosvanderVelden
SAS Super FREQ
Can you provide the code you are running (also include the xxx.sas-files)?
Is there a document from WRDS researcher that documents the steps?

Best regards, Jos
tammytt
Fluorite | Level 6

Hi Jos,

 

Thanks for your response.

 

Attached please find the two .sas files for the programs. The RollGibbsLibrary02 .sas is where the subroutine is built, and the crspGibbsBuildv01.sas is the main program.

 

This website basically describes the procedure of the programs. https://pages.stern.nyu.edu/~jhasbrou/Research/GibbsCurrent/gibbsCurrentIndex.html 

Best regards,

Tammy

JosvanderVelden
SAS Super FREQ
I suspect the error you have has to do with something in your environment. Where are you running this code? And in what order are you running the programs? Can you run both programs first RollGibbs... and attach the logfile to the track?
tammytt
Fluorite | Level 6

Hi Jos,

 

I am running the code on the WRDS(Wharton Research Data Services) cloud, since the code needs the data from wrds and that is also how the code is designed.

 

The code is run following the author's instruction (first Roll02 and then crsp01).

 

As you could see from the log file of crsp01, the error occurs at the loop on iSample 18001. I have tried to start another run from iSample 18001 - change the code to 

do iSample=18001 to nSamples;

and continues to start a new run from the new error point at some 'iSample'. I have an interesting finding: the output of the main loop 'gibbsOut' always generates 13962 rows, and a new error will occur. Do you have any idea how this problem could happen?

 

The logfile can't be attached. I will paste it in another reply.

 

Best,

Tammy

 

 

tammytt
Fluorite | Level 6

Roll02:

tammytt_0-1629221911220.png

the beginning is the setting of cloud

1          *	Gibbs estimation of the Roll model
2          
3          	IML subroutines used to estimate roll model in SAS.
4          	
5          	July 2010
6          	
7          	Joel Hasbrouck
8          _____________________________________________________________________
8        ! _______________________________
9          
10         	RollGibbs Library 02.sas
11         	
12         	Running this program compiles IML modules needed to do Gibbs
12       ! estimation of Roll model and
13         	places them in a library called IMLStor.
14         
15         	Note:
16         	There are three routines to draw (simulate) the trade direction
16       ! indicators (the q(t)).
17         	qDraw		Is written for efficiency, but it may be difficult to follow.
18         				It is a block sampler.
19         				It can only be used when c and the variance of u are fixed.
20         	qDrawSlow	Is written for clarity. It draws the q one-at-a-time,
20       ! sequentially.
21         				It can only be used when c and the variance of u are fixed.
22         	qDrawVec	This is like qDraw, but can be used when c and variance of
22       ! u are time-varying.
23         				(THIS ROUTINE HAS NOT BEEN THOROUGHLY TESTED.)
24         
25         _____________________________________________________________________
25       ! _______________________________;
26         
27         options nodate nocenter nonumber ps=70 ls=120 nomprint;
The SAS System

28         libname this '.';
NOTE: Libref THIS was successfully assigned as follows: 
      Engine:        V9 
      Physical Name: /scratch/tum/yujing_work/test_1965
29         
30         %let Infinity=1e30;
31         %let eps=1e-30;
32         
33         *______________________________________________________________________________________________
34         
35         	Define subroutines that will be used in simulations.
36          ______________________________________________________________________________________________;
37         proc iml;
NOTE: IML Ready
38         
39         start main;
40         reset storage=this.imlstor;
41         store module=_all_;
42         remove module=main;
43         show storage;
44         finish main;
NOTE: Module MAIN defined.
45         
46         
47         *____________________________________________________________________________________________________
48         
49         	RollGibbsBeta: Estimate Roll model with Gibbs sampler
50         	
51         	The return argument is parmOut[nSweeps,3]
52         	Column 1:	c
53         	Column 2: beta
54         	Column 3: varu
55         ____________________________________________________________________________________________________;
56         
57         
58         %let Infinity=1e30;
59         %let eps=1e-30;
60         
61         start RollGibbsBeta(parmOut, p, pm, q, nSweeps, regDraw, varuDraw, qDraw, varuStart, cStart, betaStart,
61       ! printLevel);
62         	
62       !  nObs = nrow(p);
63         	
63       !  if nrow(q)^=nObs | nrow(pm)^=nObs then do;
64         		
64       !   print 'RollGibbsBeta length mismatch';
65         		
65       !   return;
66         	
66       !  end;
67         	
67       !  dp = p[2:nObs] - p[1:(nObs-1)];
68         	
69         	
69       !  if qDraw then do;
69       !                    	*	Initialize qs to price sign changes;
70         		
70       !   qInitial = {1} // sign(dp);
71         		
71       !   qInitial = qInitial # (q^=0);
71       !                                	*	Only initialize nonzero elements of q;
72         		
72       !   q = qInitial;
73         	
73       !  end;
74         
75         	
75       !  if varuStart<=0 then varuStart = 0.001;
The SAS System

76         	
76       !  varu = varuStart;
77         	
77       !  if cStart<=0 then cStart=0.01;
78         	
78       !  c = cStart;
79         	
79       !  if betaStart<=0 then betaStart=1;
80         	
80       !  beta = betaStart;
81         	
82         	
82       !  parmOut = j(nSweeps,3,.);
83         
84         	
84       !  do sweep=1 to nSweeps;
85         
86         		
86       !   dq =  q[2:nObs] - q[1:(nObs-1)];
87         		
87       !   dpm = pm[2:nObs] - pm[1:(nObs-1)];
88         
89         		
89       !   if regDraw then do;
90         			
90       !    priorMu = {0,1};
91         			
91       !    postMu = priorMu;
92         			
92       !    priorCov = diag({1,2});
93         			
93       !    postCov = priorCov;
94         			
94       !    X = colvec(dq) || colvec(dpm);
95         			
95       !    rc = BayesRegressionUpdate(priorMu, priorCov, dp, X, varu, postMu, postCov);
96         			
96       !    if printLevel>=2 then print postMu postCov;
97         			
97       !    coeffLower={0,-&Infinity};
98         			
98       !    coeffUpper=j(2,1,&Infinity);
99         			
99       !    coeffDraw = mvnrndT(postMu, postCov, coeffLower, coeffUpper);
99       !                                                                  /* */
100        			
100      !    if printLevel>=2 then print coeffDraw;
101        			
101      !    c = coeffDraw[1];
102        			
102      !    beta = coeffDraw[2];
103        		
103      !   end;
104        
105        		
105      !   if varuDraw then do;
106        			
106      !    u = dp - c*dq - beta*dpm;
107        			
107      !    priorAlpha = 1.e-12;
108        			
108      !    priorBeta = 1.e-12;
109        			
109      !    postAlpha = .;
110        			
110      !    postBeta = .;
111        			
111      !    rc = BayesVarianceUpdate(priorAlpha, priorBeta, u, postAlpha, postBeta);
The SAS System

112        			
112      !    x = (1/postBeta) * rand('gamma',postAlpha);
113        			
113      !    varu = 1/x;
114        			
114      !    sdu = sqrt(varu);
115        			
115      !    if printLevel>=2 then print varu;
116        		
116      !   end;
117        
118        		
118      !   if qDraw then do;
119        			
119      !    qDrawPrintLevel = 0;
120        			
120      !    p2 = p - beta*pm;
121        			
121      !    call qDraw(p2, q, c, varu, qDrawPrintLevel);
122        		
122      !   end;
123        		
124        		
124      !   parmOut[sweep, 1]=c;
125        		
125      !   parmOut[sweep, 2] = beta;
126        		
126      !   parmOut[sweep, 3] = varu;
127        		
128        end;
129        finish RollGibbsBeta;
NOTE: Module ROLLGIBBSBETA defined.
130        
131        
132        *______________________________________________________________________________________________
133        
134        call qDraw(p, q, c, varu, printLevel)	makes new draws for q
135        
136        parameters:
137        p		column vector of trade prices
138        q		column vector of qs (REPLACED ON RETURN)
139        c		cost parameter
140        varu	variance of disturbance
141        printLevel	used to generate debugging output
142        
143        _______________________________________________________________________________________________;
144        start qDraw(p, q, c, varu, printLevel);
145        if nrow(p)^=nrow(q) then do;
146        	
146      !  print "qDraw. p and q are of different lengths. p is " nrow(p) "x" ncol(p) "; q is " nrow(q) "x" ncol(q);
147        	
147      !  abort;
148        end;
149        if nrow(c)^=1 then do;
150        	
150      !  print "qDraw. c should be a scalar. It is " nrow(c) "x" ncol(c);
151        	
151      !  abort;
152        end;
153        if nrow(varu)^=1 then do;
154        	
154      !  print "qDraw. varu should be a scalar. It is " nrow(varu) "x" ncol(varu);
155        	
155      !  abort;
156        end;
157        if printlevel>0 then print p q;
158        qNonzero = q^=0;
159        q = q || q;
The SAS System

160        p2 = p || p;
161        nSkip = 2;
162        modp = colvec( mod(1:nrow(p),nSkip) );
163        ru = uniform(j(nrow(p),1,0));
164        do iStart=0 to (nSkip-1);
165        	
165      !  if printLevel>0 then print iStart [l="" r="qDraw. iStart:" f=1.];
166        	
166      !  k = modp=iStart;
167        	
167      !  jnz = loc( k & qNonzero );
167      !                            	*	These are the q's we'll be drawing.;
168        	
168      !  if printLevel>0 then print jnz [l="Drawing q's for t=" f=3.];
169        	
169      !  q[jnz,1] = 1;
170        	
170      !  q[jnz,2] = -1;
171        	
171      !  cq = c*q;
172        	
172      !  v = p2 - cq;
173        	
173      !  u = v[2:nrow(v),] - v[1:(nrow(v)-1),];
174        	
174      !  if printLevel>0 then print u [l="u:"];
175        	
175      !  s=(u##2)/(2*varu);
176        	
176      !  if printLevel>0 then print s [l="s"];
177        	
177      !  sSum = (s//{0 0}) + ({0 0}//s);
178        	
178      !  if printLevel>0 then print sSum [l="sSum (before reduction)"];
179        	
179      !  sSum = sSum[jnz,];
180        	
180      !  if printLevel>0 then print sSum [l="sSum (after reduction)"];
181        	
181      !  logOdds = sSum[,2] - sSum[,1];
182        	*	Make sure that we won't get overflow when we call exp();
183        	
183      !  logOkay = logOdds<500;
184        	
184      !  Odds = exp(logOkay#logOdds);
185        	
185      !  pBuy = Odds/(1+Odds);
186        	
186      !  pBuy = logOkay#pBuy + ^logOkay;
187        	
187      !  if printLevel>0 then print pBuy [f=e10.];
188        	
188      !  qknz = 1 - 2*(ru[jnz]>pBuy);
189        	
189      !  q[jnz,1] = qknz;
190        	
190      !  if istart<(nSkip-1) then q[jnz,2] = qknz;
191        end;
192        q = q[,1];
193        finish qDraw;
NOTE: Module QDRAW defined.
194        
195        
196        *______________________________________________________________________________________________
197        
198        call qDrawVec(p, q, c, varu, printLevel)	makes new draws for q
199        
200        parameters:
The SAS System

201        p		column vector of trade prices
202        q		column vector of qs (REPLACED ON RETURN)
203        c		cost parameter (either a scalar or a Tx1 vector)
204        varu	variance of disturbance (either a scalar or a T-1 x 1 vector)
205        printLevel	used to generate debugging output
206        
207        _______________________________________________________________________________________________;
208        start qDrawVec(p, q, c, varu, printLevel);
209        if nrow(p)^=nrow(q) then do;
210        	
210      !  print "qDraw. p and q are of different lengths. p is " nrow(p) "x" ncol(p) "; q is " nrow(q) "x" ncol(q);
211        	
211      !  abort;
212        end;
213        if nrow(c)^=1 & nrow(c)^=nrow(p) then do;
214        	
214      !  print "qDraw. p is  " nrow(p) "x" ncol(p) "; c is " nrow(c) "x" ncol(c);
215        	
215      !  abort;
216        end;
217        if nrow(varu)^=1 & nrow(varu)^=(nrow(p)-1) then do;
218        	
218      !  print "qDraw. p is  " nrow(p) "x" ncol(p) "; varu is " nrow(varu) "x" ncol(varu) "(should  be T-1 x 1)";
219        	
219      !  abort;
220        end;
221        if printlevel>0 then print p q;
222        qNonzero = q^=0;
223        if nrow(c)^=1 then c2 = c || c;
224        q = q || q;
225        p2 = p || p;
226        if nrow(varu)^=1 then varu2 = 2*(varu || varu);
227        else varu2 = 2*varu;
228        nSkip = 2;
229        modp = colvec( mod(1:nrow(p),nSkip) );
230        do iStart=0 to (nSkip-1);
231        	
231      !  if printLevel>0 then print iStart [l="" r="qDraw. iStart:" f=1.];
232        	
232      !  k = modp=iStart;
233        	
233      !  jnz = loc( k & qNonzero );
233      !                            	*	These are the q's we'll be drawing.;
234        	
234      !  if printLevel>0 then print jnz [l="Drawing q's for t=" f=3.];
235        	
235      !  q[jnz,1] = 1;
236        	
236      !  q[jnz,2] = -1;
237        	
237      !  if nrow(c)=1 then cq = c*q;
238        	
238      !  else cq = c2#q;
239        	
239      !  v = p2 - cq;
240        	
240      !  u = v[2:nrow(v),] - v[1:(nrow(v)-1),];
241        	
241      !  if printLevel>0 then print u [l="u:"];
242        	
242      !  if nrow(varu)=1 then s=(u##2)/(2*varu);
243        	
243      !  else s = (u##2)/varu2;
244        	
244      !  if printLevel>0 then print s [l="s"];
245        	
245      !  sSum = (s//{0 0}) + ({0 0}//s);
246        	
The SAS System

246      !  if printLevel>0 then print sSum [l="sSum (before reduction)"];
247        	
247      !  sSum = sSum[jnz,];
248        	
248      !  if printLevel>0 then print sSum [l="sSum (after reduction)"];
249        	
249      !  logOdds = sSum[,2] - sSum[,1];
250        	*	Make sure that we won't get overflow when we call exp();
251        	
251      !  logOkay = logOdds<500;
252        	
252      !  Odds = exp(logOkay#logOdds);
253        	
253      !  pBuy = Odds/(1+Odds);
254        	
254      !  pBuy = logOkay#pBuy + ^logOkay;
255        	
255      !  if printLevel>0 then print pBuy [f=e10.];
256        	
256      !  ru = uniform( j(nrow(pBuy),1) );
257        	
257      !  qknz = 1 - 2*(ru>pBuy);
258        	
258      !  q[jnz,1] = qknz;
259        	
259      !  if istart<(nSkip-1) then q[jnz,2] = qknz;
260        end;
261        q = q[,1];
262        finish qDrawVec;
NOTE: Module QDRAWVEC defined.
263        
264        *______________________________________________________________________________________________
265        
266        call qDrawSlow(p, q, c, varu)	makes new draws for q
267        inputs:
268        p		column vector of trade prices
269        q		column vector of qs (REPLACED ON RETURN)
270        c		cost parameter
271        varu	variance of disturbance
272        
273         ______________________________________________________________________________________________;
274        start qDrawSlow(p, q, c, varu, PrintLevel);
275        T = nrow(p);
276        reset noname spaces=1;
277        if PrintLevel>0 then do;
278        	
278      !  print "qDraw" T [r="T:"] varu [l="" r="varu:"] c [r="c:"];
279        	
279      !  print (t(q));
280        end;
281        do s=1 to T;
282        	
282      !  if PrintLevel>=2 then print s [r="s:"];
283        	
283      !  if q[s]=0 then goto sNext;
283      !                            	* Don't make a draw if q=0 (p is a quote midpoint);
284        	
284      !  prExp = j(1,2,0);
285        	
285      !  if s<T then do;
286        		
286      !   uAhead = (p[s+1] - p[s]) + c*{1 -1} - c*q[s+1];
287        		
287      !   prExp = prExp - (uAhead##2)/(2*varu);
288        		
288      !   if PrintLevel>=2 then print s [format=5. r='s'] (p[s])   [r='p[s]  '] (p[s+1]) [r='p[s+1]'] uAhead
288      ! [r='uAhead'];
289        	
The SAS System

289      !  end;
290        	
290      !  if s>1 then do;
291        		
291      !   uBack = (p[s] -  p[s-1]) + c*q[s-1] - c*{1 -1};
292        		
292      !   prExp = prExp - (uBack##2)/(2*varu);
293        		
293      !   if PrintLevel>=2 then print s [format=5. r='s'] (p[s-1]) [r='p[s-1]'] (p[s])   [r='p[s]  '] uBack
293      ! [r='uBack '];
294        	
294      !  end;
295        	
295      !  logOdds = prExp[1]-prExp[2];
296        	
296      !  if PrintLevel>=2 then print logOdds [r='logOdds (in favor of buy)'];
297        	
297      !  if abs(logOdds)>100 then q[s]=sign(logOdds);
298        	
298      !  else do;
299        		
299      !   pBuy = 1 - 1/(1+exp(logOdds));
300        		
300      !   q[s] = 1 - 2*(rand('uniform')>PBuy);
301        	
301      !  end;
302        sNext:
302      !        end;
303        reset name;
304        return;
305        finish qDrawSlow;
NOTE: Module QDRAWSLOW defined.
306        
307        
308        
309        *______________________________________________________________________________________________
310        
311        rc = BayesVarianceUpdate(priorAlpha, priorBeta, u, postAlpha, postBeta)
312        	updates the variance posterior (inverted gamma)
313        inputs:
314        priorAlpha and priorBeta
315        u	vector of estimated disturbances
316        postAlpha and postBeta are updated on return to the posterior values
317        _______________________________________________________________________________________________;
318        start BayesVarianceUpdate(priorAlpha, priorBeta, u, postAlpha, postBeta);
319        postAlpha = priorAlpha + nrow(u)/2;
320        postBeta = priorBeta + (u##2)[+]/2;
321        return (0);
322        finish BayesVarianceUpdate;
NOTE: Module BAYESVARIANCEUPDATE defined.
323        
324        *______________________________________________________________________________________________
325        
326        rc = BayesRegressionUpdate(priorMu, priorCov, y, X, dVar, postMu, postCov)
327        	computes coefficient posteriors for normal Bayesian regression model
328        inputs:
329        priorMu and priorCov	coefficient priors (mean and covariance)
330        y	column vector of l.h.s. values
331        X	matrix of r.h.s. variables
332        dVar	error variance (taken as given)
333        postMu and postCov 	coefficient posteriors.
334        rc is a return code (0 if okay)
335        _______________________________________________________________________________________________;
336        start BayesRegressionUpdate(priorMu, priorCov, y, X, dVar, postMu, postCov);
337        if ncol(priorMu)^=1 then do;
338        	
338      !  print "BayesRegressionUpdate. priorMu is " nrow(priorMu) "x" ncol(priorMu) " (should be a column vector)";
339        	
The SAS System

339      !  return(-1);
340        end;
341        if nrow(X)<ncol(X) then do;
342        	
342      !  print "BayesRegressionUpdate. X is " nrow(X) "x" ncol(X);
343        	
343      !  return (-1);
344        end;
345        if nrow(X)^=nrow(y) | ncol(y)^=1 then do;
346        	
346      !  print "BayesRegressionUpdate. X is " nrow(X) "x" ncol(X) "; y is " nrow(y) "x" ncol(y);
347        	
347      !  return (-1);
348        end;
349        if nrow(priorMu)^=ncol(X) then do;
350        	
350      !  print "BayesRegressionUpdate. X is " nrow(X) "x" ncol(X) "; priorMu is " nrow(priorMu) "x" ncol(priorMu)
351        		" (not conformable)";
352        	
352      !  return (-1);
353        end;
354        if nrow(priorCov)^=ncol(priorCov) | nrow(priorCov)^=nrow(priorMu) then do;
355        	
355      !  print "BayesRegressionUpdate. priorMu is " nrow(X) "x" ncol(X) "; priorCov is " nrow(priorCov) "x"
355      ! ncol(priorCov);
356        	
356      !  return (-1);
357        end;
358        
359        covi = inv(priorCov);
360        Di = (1/dVar)*(t(X)*X) + covi;
361        D = inv(Di);
362        dd = (1/dVar)*t(X)*y + covi*priorMu;
363        postMu = D*dd;
364        postCov = D;
365        return (0);
366        finish BayesRegressionUpdate;
NOTE: Module BAYESREGRESSIONUPDATE defined.
367        
368        *______________________________________________________________________________________________
369        
370        RandStdNormT(zlow,zhigh) returns a random draw from the standard normal distribution
371        truncated to the range (zlow, zhigh).
372        _______________________________________________________________________________________________;
373        start RandStdNormT(zlow,zhigh);
374        if zlow=-&Infinity & zhigh=&Infinity then return (normal(seed));
375        PROBNLIMIT = 6;
376        if zlow>PROBNLIMIT & (zhigh=&Infinity | zhigh>PROBNLIMIT) then return (zlow+100*&eps);
377        if zhigh<-PROBNLIMIT & (zlow=-&Infinity | zlow<-PROBNLIMIT) then return (zhigh-100*&eps);
378        if zlow=-&Infinity then plow=0;
379        else plow = probnorm(zlow);
380        if zhigh=&Infinity then phigh=1;
381        else phigh = probnorm(zhigh);
382        p = plow + rand('uniform')*(phigh-plow);
382      !                                           /* should be the source of error */
383        if p=1 then return (zlow + 100*eps);
384        if p=0 then return (zhigh - 100*eps);
385        return (probit(p));
386        finish RandStdNormT;
NOTE: Module RANDSTDNORMT defined.
387        
388        *______________________________________________________________________________________________
389        
390        mvnrndT(mu, cov, vLower, vUpper) returns a random draw (a vector) from a multivariate normal
391        distribution with mean mu and covariance matrix cov, truncated to the range (vLower, vUpper).
392        vLower and vUpper are vectors (conformable with mu) that specify the upper and lower truncation
393        points for each component.
394        _______________________________________________________________________________________________;
The SAS System

395        start mvnrndT(mu, cov, vLower, vUpper);
396        f = t(root(cov));
397        n = nrow(mu)*ncol(mu);
398        eta = j(n,1,0);
399        low = (vLower[1]-mu[1])/f[1,1];
400        high = (vUpper[1]-mu[1])/f[1,1];
401        eta[1] = RandStdNormT(low,high);
402        do k=2 to n;
403        	
403      !  etasum = f[k,1:(k-1)]*eta[1:(k-1)];
404        	
404      !  low = (vLower[k]-mu[k]-etasum)/f[k,k];
405        	
405      !  high = (vUpper[k]-mu[k]-etasum)/f[k,k];
406        	
406      !  eta[k] = RandStdNormT(low,high);
407        end;
408        return (colvec(mu)+f*eta);
409        finish mvnrndT;
NOTE: Module MVNRNDT defined.
410        
411        run;
NOTE: New storage library = THIS.IMLSTOR
412        quit;
NOTE: Exiting IML.
NOTE: Storage library THIS.IMLSTOR closed.
NOTE: The PROCEDURE IML printed page 1.
NOTE: PROCEDURE IML used (Total process time):
      real time           0.08 seconds
      cpu time            0.04 seconds
      

413        
414        
415        *______________________________________________________________________________________________
416        
417        	Test the routines
418         ______________________________________________________________________________________________;
419        proc iml;
NOTE: IML Ready
420        	
420      !  start main;
421        	
421      !  call streaminit(1234);
422        	
422      !  reset storage=this.imlstor;
423        	
423      !  load;
424        
425        	
425      !  print 'Test of RandStdNorm:';
426        	
426      !  x =RandStdNormT(.5,1);
427        	
427      !  print x;
428        	
429        	
429      !  print 'Test of mvnrndT:';
430        	
430      !  mu={1 2};
431        	
431      !  cov={1 .5, .5 1};
432        	
432      !  vLower = {0 0};
433        	
433      !  vUpper = {&infinity &infinity};
434        	
434      !  x = mvnrndT(mu, cov, vLower, vUpper);
The SAS System

435        	
435      !  print x;
436        	
437        	
437      !  print 'Test of Bayesian normal regression update';
438        	
438      !  nObs = 1000;
439        	
439      !  nv = 2;
440        	
440      !  x = j(nObs,nv);
441        	
441      !  u = j(nObs,1);
442        	
442      !  sdu = 2;
443        	
443      !  do t=1 to nObs;
444        		
444      !   x[t,1] = 1;
445        		
445      !   u = sdu*rand('normal');
446        		
446      !   do i=2 to nv;
447        			
447      !    x[t,i] = rand('normal');
448        		
448      !   end;
449        	
449      !  end;
450        	
450      !  y = x * colvec((1+(1:nv))) + u;
451        	
451      !  priorMu={0,0};
452        	
452      !  priorCov={&Infinity 0, 0 &Infinity};
453        	
453      !  dVar = sdu*sdu;
454        	
454      !  postMu = 0;
455        	
455      !  postCov = 0;
456        	
456      !  rc = BayesRegressionUpdate(priorMu, priorCov, y, X, dVar, postMu, postCov);
457        	
457      !  print postMu;
458        	
458      !  print postCov;
459        	
459      !  priorAlpha = 1.e-6;
460        	
460      !  priorBeta = 1.e-6;
461        	
461      !  postAlpha = 0;
462        	
462      !  postBeta = 0;
463        	
463      !  u = y - X*postMu;
464        	
464      !  rc = BayesVarianceUpdate(priorAlpha, priorBeta, u, postAlpha, postBeta);
465        	
465      !  print postAlpha;
466        	
466      !  print postBeta;
467        	
467      !  finish main;
NOTE: Module MAIN defined.
468        	
468      !  run;
The SAS System

NOTE: New storage library = THIS.IMLSTOR
NOTE: Module COLVEC loaded from the storage SASHELP.IMLMLIB.
469        	quit;
NOTE: Exiting IML.
NOTE: Storage library THIS.IMLSTOR closed.
NOTE: The PROCEDURE IML printed page 2.
NOTE: PROCEDURE IML used (Total process time):
      real time           0.02 seconds
      cpu time            0.01 seconds
      

470        	
471        proc iml;
NOTE: IML Ready
472        	
472      !  start main;
473        	
473      !  call streaminit(1234);
474        	
474      !  reset storage=this.imlstor;
475        	
475      !  load;
476        
477        	
477      !  nObs = 250;
477      !              /* the program WORKS if NOBS>=5, but fails for small NOBS */
478        	
478      !  sdm = sqrt(0.3##2 / 250);
479        	
479      !  sdu = sqrt(0.2##2 / 250);
480        	
480      !  print sdu;
481        	
481      !  u = j(nObs,1,.);
482        	
482      !  v = j(nObs,1,.);
483        	
483      !  call randseed(12345);
484        	
484      !  call randgen(u,'normal',0,sdu);
485        	
485      !  call randgen(v,'normal',0,sdm);
486        	
486      !  beta = 1.1;
487        	
487      !  c = .04;
488        	
488      !  q = j(nObs,1,.);
489        	
489      !  call randgen(q,'uniform');
490        	
490      !  q = sign(q-0.5);
491        	
491      !  probZero = 0.6;
492        	
492      !  z = j(nObs,1,.);
493        	
493      !  call randgen(z,'uniform');
494        	
494      !  q = (z>probZero)#q;
495        	
495      !  p = j(nObs,1,.);
496        	
496      !  m = 0;
497        	
497      !  do t=1 to nObs;
498        		
498      !   m = m + beta*v[t] + u[t];
The SAS System

499        		
499      !   p[t] = m + c*q[t];
500        	
500      !  end;
501        	
501      !  pM = t(cusum(t(v)));
502        	
503        	
503      !  nSweeps = 1000;
504        	
504      !  regDraw = 1;
505        	
505      !  varuDraw = 1;
506        	
506      !  qDraw = 1;
507        	
507      !  nDrop = 200;
508        	
509        	
509      !  call RollGibbsBeta(parmOut, p, pm, q, nSweeps, regDraw, varuDraw, qDraw, 0,0,0,0);
510        	
511        	
511      !  p2 = parmOut[(nDrop+1):nSweeps,];
512        	
512      !  p2 = p2 || sqrt(p2[,3]);
513        	
513      !  pm = p2[+,]/(nSweeps-nDrop);
514        	
514      !  print pm;
515        
516        
517        	
517      !  finish main;
NOTE: Module MAIN defined.
518        run;
NOTE: New storage library = THIS.IMLSTOR
NOTE: Module COLVEC loaded from the storage SASHELP.IMLMLIB.
519        quit;
NOTE: Exiting IML.
NOTE: Storage library THIS.IMLSTOR closed.
NOTE: The PROCEDURE IML printed page 3.
NOTE: PROCEDURE IML used (Total process time):
      real time           0.16 seconds
      cpu time            0.16 seconds
      

519      !      		
520        	

NOTE: SAS Institute Inc., SAS Campus Drive, Cary, NC USA 27513-2414
NOTE: The SAS System used:
      real time           1.08 seconds
      cpu time            0.43 seconds
      
tammytt
Fluorite | Level 6

The log file is longer than I thought, I have uploaded it in google drive. Sorry for the inconvenience.

https://drive.google.com/drive/folders/14yrCDn2cJlr47WvKcuA9D4ArrlarFVv3?usp=sharing 

sas-innovate-white.png

Register Today!

Join us for SAS Innovate 2025, our biggest and most exciting global event of the year, in Orlando, FL, from May 6-9.

 

Early bird rate extended! Save $200 when you sign up by March 31.

Register now!

How to Concatenate Values

Learn how use the CAT functions in SAS to join values from multiple variables into a single value.

Find more tutorials on the SAS Users YouTube channel.

SAS Training: Just a Click Away

 Ready to level-up your skills? Choose your own adventure.

Browse our catalog!

Discussion stats
  • 6 replies
  • 1466 views
  • 0 likes
  • 2 in conversation