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;
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
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
Roll02:
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
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
It's finally time to hack! Remember to visit the SAS Hacker's Hub regularly for news and updates.
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.
Ready to level-up your skills? Choose your own adventure.