Hi Koen,
Thanks for your help!
I haven't run it successfully before. But I do get the codes of testing the subroutines, which works just fine as following.
*______________________________________________________________________________________________
Test the routines
______________________________________________________________________________________________;
proc iml;
start main;
call streaminit(1234);
reset storage=this.imlstor;
load;
print 'Test of RandStdNorm:';
x =RandStdNormT(.5,1);
print x;
print 'Test of mvnrndT:';
mu={1 2};
cov={1 .5, .5 1};
vLower = {0 0};
vUpper = {&infinity &infinity};
x = mvnrndT(mu, cov, vLower, vUpper);
print x;
print 'Test of Bayesian normal regression update';
nObs = 1000;
nv = 2;
x = j(nObs,nv);
u = j(nObs,1);
sdu = 2;
do t=1 to nObs;
x[t,1] = 1;
u = sdu*rand('normal');
do i=2 to nv;
x[t,i] = rand('normal');
end;
end;
y = x * colvec((1+(1:nv))) + u;
priorMu={0,0};
priorCov={&Infinity 0, 0 &Infinity};
dVar = sdu*sdu;
postMu = 0;
postCov = 0;
rc = BayesRegressionUpdate(priorMu, priorCov, y, X, dVar, postMu, postCov);
print postMu;
print postCov;
priorAlpha = 1.e-6;
priorBeta = 1.e-6;
postAlpha = 0;
postBeta = 0;
u = y - X*postMu;
rc = BayesVarianceUpdate(priorAlpha, priorBeta, u, postAlpha, postBeta);
print postAlpha;
print postBeta;
finish main;
run;
quit;
proc iml;
start main;
call streaminit(1234);
reset storage=this.imlstor;
load;
nObs = 250;
sdm = sqrt(0.3##2 / 250);
sdu = sqrt(0.2##2 / 250);
print sdu;
u = j(nObs,1,.);
v = j(nObs,1,.);
call randseed(12345);
call randgen(u,'normal',0,sdu);
call randgen(v,'normal',0,sdm);
beta = 1.1;
c = .04;
q = j(nObs,1,.);
call randgen(q,'uniform');
q = sign(q-0.5);
probZero = 0.6;
z = j(nObs,1,.);
call randgen(z,'uniform');
q = (z>probZero)#q;
p = j(nObs,1,.);
m = 0;
do t=1 to nObs;
m = m + beta*v[t] + u[t];
p[t] = m + c*q[t];
end;
pM = t(cusum(t(v)));
nSweeps = 1000;
regDraw = 1;
varuDraw = 1;
qDraw = 1;
nDrop = 200;
call RollGibbsBeta(parmOut, p, pm, q, nSweeps, regDraw, varuDraw, qDraw, 0,0,0,0);
p2 = parmOut[(nDrop+1):nSweeps,];
p2 = p2 || sqrt(p2[,3]);
pm = p2[+,]/(nSweeps-nDrop);
print pm;
finish main;
run;
quit; 
Thanks so much for your help!
Another information is although the module RollGibbsBeta didn't run till the end successfully, it did produce (a part of) "gibbsOut" file with 13962 rows and 7 columns, and the results seems quite appropriate with the output in c, beta, varu and sdu (at least looks like so). Could that mean the module RollGibbsBeta works fine at the beginning and some inappropriate data in the middle make it fail?
the gibbsOut output:
the main program:
*____________________________________________________________________________________________________
This proc loops over all samples, calling the estimation routines.
_____________________________________________________________________________________________________;
options nonotes;
proc iml;
start main;
call streaminit(1234);
reset storage=this.imlstor; * This contains the subroutines built in RollGibbsLibrary02.sas;
load;
reset printadv=1 log;
use temp.pys;
read all var {permno year kSample} where(nTradeDays>=60) into sample [colname=colSample]; *my correction;
nSamples = nrow(sample);
print 'nSamples=' nSamples;
permno = sample[,1];
year = sample[,2];
kSample = sample[,3];
outSet = j(1,7,.);
varnames ={'permno','year','kSample','c','beta','varu','sdu'};
create this.gibbsOut from outSet [colname=varnames];
do iSample=1 to nSamples;
thisPermno = permno[iSample];
thisYear = year[iSample];
thisKSample = kSample[iSample];
if mod(iSample,500)=1 then do;
t = time();
ctime = putn(t,'time.');
print ctime iSample '/' nSamples ': ' thisPermno thisYear thisKSample;
end;
if thisPermno>=&startPermno & thisYear>=1990 then do;
use temp.dsf where (permno=thisPermno & year=thisYear & kSample=thisKSample);
read all var {p pm q} into x [colname=colx];
nSweeps = 1000;
regDraw = 1;
varuDraw = 1;
qDraw = 1;
nDrop = 200;
call RollGibbsBeta(parmOut, x[,1],x[,2],x[,3], nSweeps, regDraw, varuDraw, qDraw, 0,0,0,0);
p2 = parmOut[(nDrop+1):nSweeps,];
p2 = p2 || sqrt(p2[,3]);
pm = p2[+,]/(nSweeps-nDrop);
outset = thisPermno || thisYear || thisKSample || pm;
*print outset;
setout this.gibbsOut;
append from outSet;
end;
end;
finish main;
run;
quit;
options notes;
proc print data=this.gibbsOut (obs=50);
run;
data this.crspGibbs2009;
set this.gibbsOut;
run;
Multiple Linear Regression in SAS
Learn how to run multiple linear regression models with and without interactions, presented by SAS user Alex Chaplin.
Find more tutorials on the SAS Users YouTube channel.
-
2628 views
-
2 likes
-
-
