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
... View more