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