Random numbers seed setting

Hi,

I am using randn('seed', 1) at the beginning of my code for a simulation with 1000 replications. This is to make sure that if I rerun the simulation I'll get same results every time. But I am getting slightly different results in different runs. Can anyone help? I am using Matlab 2012a.

Answer 1

[Azzi Abdelmalek]

Use

a=rng;
out=randn
rng(a)
out=randn

[Rabeya]

Thanks. Do you suspect the randn('seed',1) is not working in Matlab 2012a?

[Azzi Abdelmalek]

It's working for me

randn('seed',1)
randn
randn('seed',1)
randn

[Rabeya]

Thanks again. Actually I am facing a very strange problem. Sometimes it works, sometimes it does not, for me. Its fine with a random generation, but when I run 1000 replications and report the mean, median of the results, they are different in different runs. Do you have any idea what is going on there?

[Azzi Abdelmalek]

can you post the part of your code using randn

[Rabeya]

clear
clc
randn('seed',1);
rep =1000; % number of replications
k=2; % number of parameters
onesr = ones(rep,1);
beta_true=[0;0.05];

% for data:
obsvec=100;%[5;10;50;100];%;1000];
obsval=length(obsvec);
expand=26; % each cohort observed over 26 years

Cvec=30;%[27;30;40;60];%[60;12;6]; % no of cohorts
Gvec=15;%[3;10;30];%[3;3;4;6];%[4;10;40]; %[6;2]; % no of groups
Cval=length(Cvec);
Gval = length(Gvec);

ci = 0;
while ci < Cval %takes different number of cohorts
  ci = ci+1;
  C = Cvec(ci);
  %G = Gvec(ti);
  oi=0;
  while oi<obsval
      oi=oi+1;
      obs=obsvec(oi);
      gi=0;
      while gi < Gval %takes different number of groups
          gi=gi+1;
          G=Gvec(gi);
          %  storage for the beta-estimator values
          beta2sls=zeros(rep,k);
          se_beta2sls=zeros(rep,k);

            betaGMM=zeros(rep,k); 
            se_betaGMM=zeros(rep,k);
            J_stat_GMM=zeros(rep,1);
            pval_GMM=zeros(rep,1);

            beta0=zeros(rep,1);
            betaEL=zeros(rep,k);
            se_betaEL=zeros(rep,k);
            bias_crrctdEL=zeros(rep,k);
            se_bias_crrctdEL=zeros(rep,k);
            lik_EL=zeros(rep,1);
            J_statEL=zeros(rep,1);
            pvalEL=zeros(rep,1);

            betaET=zeros(rep,k);
            se_betaET=zeros(rep,k);
            bias_crrctdET=zeros(rep,k);
            se_bias_crrctdET=zeros(rep,k);
            lik_ET=zeros(rep,1);
            J_statET=zeros(rep,1);
            pvalET=zeros(rep,1);

            c_GMM=zeros(rep,1);
            c_EL=zeros(rep,1);
            c_ET=zeros(rep,1);

            wald_EL=zeros(rep,1);
            pval_wald_EL=zeros(rep,1);
            c_wald_EL=zeros(rep,1);
            LM_EL=zeros(rep,1);
            pval_LM_EL=zeros(rep,1);
            c_LM_EL=zeros(rep,1);

            wald_ET=zeros(rep,1);
            pval_wald_ET=zeros(rep,1);
            c_wald_ET=zeros(rep,1);
            LM_ET=zeros(rep,1);
            pval_LM_ET=zeros(rep,1);
            c_LM_ET=zeros(rep,1);

            bar_ug=zeros(G,rep);
            sigma2_ug=zeros(G,rep);
            ug_3=zeros(G,rep);
            xu_g=zeros(G,k-1,rep);
            N=zeros(rep,1);
            se=[2;1.8;1.75;1.71;1.6;1.5;1.45;1.41;1.4;1.36;1.2;1.11;1.05;1.05;1.02;1;0.99;0.95;0.91;0.9;0.86;0.88;0.84;0.8;0.84;0.7;0.66;0.5;0.43;0.39];

            r = 0;
            while r < rep % replication
                r = r+1;
                data=Data_RCSabilhet_step(C,obs,expand,se);
                %data=pseudo2_data(T,obs,expand); % calling the data
                [x,y,u,group,ng]=xyu_RCSabilhet_G15_withcons(data); % forming the variables
                N(r,1)=sum(ng);
                [bar_ug(:,r), sigma2_ug(:,r), ug_3(:,r),xu_g(:,:,r)]=moment_u(x,u,ng,group);
                z=dummyvar(group); % group dummies, needed for GMM estimation

                % GMM estimation
                [beta2sls(r,:),se_beta2sls(r,:), iter, betaGMM(r,:),se_betaGMM(r,:), J_stat_GMM(r,:), pval_GMM(r,:)]=igmm_grouped_data(y,x,z,ng);
                if  pval_GMM(r,:)>=0.05 % at 5% sig level
                    c_GMM(r,:)=1; % number of rejections
                end    

                % EL estimators:
                [betaEL(r,:),lik_EL(r,:),lambda,p]=EL_optim_owen(x,y,G,ng,betaGMM(r,:)');
                se_betaEL(r,:)=diag(se_beta_GEL(betaEL(r,:)',x,y,group,ng,p));

                bias_EL=bias_GEL(x,y,G,ng,betaEL(r,:)',-2,p);
                bias_crrctdEL(r,:)=betaEL(r,:)-bias_EL';
                se_bias_crrctdEL(r,:)=diag(se_beta_GEL(bias_crrctdEL(r,:)',x,y,group,ng,p));
                J_statEL(r,:)=2*lik_EL(r,:); % overidentification test stat 
                pvalEL(r,:)=1-chi2cdf(J_statEL(r,:),G-k); % pvalue of the test statistic
                if  pvalEL(r,:)>=0.05 % at 5% sig level
                    c_EL(r,:)=1; % number of rejections
                end
                wald_EL(r,:)=wald_test(x,y,ng,G,betaEL(r,:)',p);
                pval_wald_EL(r,:)=1-chi2cdf(wald_EL(r,:),G-k);
                if  pval_wald_EL(r,:)>=0.05 % at 5% sig level
                    c_wald_EL(r,:)=1; % number of rejections
                end
                LM_EL(r,:)=LM_test(x,y,ng,G,betaEL(r,:)',lambda,p);
                pval_LM_EL(r,:)=1-chi2cdf(LM_EL(r,:),G-k);
                if  pval_LM_EL(r,:)>=0.05 % at 5% sig level
                    c_LM_EL(r,:)=1; % number of rejections
                end
                % ET estimators:
                [betaET(r,:),lik_ET(r,:),lambda,p_ET]=ET_optim(x,y,G,ng,betaGMM(r,:)');
                se_betaET(r,:)=diag(se_beta_GEL(betaET(r,:)',x,y,group,ng,p_ET));

                bias_ET=bias_GEL(x,y,G,ng,betaET(r,:)',-1,p_ET);
                bias_crrctdET(r,:)=betaET(r,:)-bias_ET';
                se_bias_crrctdET(r,:)=diag(se_beta_GEL(bias_crrctdET(r,:)',x,y,group,ng,p_ET));
                J_statET(r,:)=(2*N(r,1)+2*lik_ET(r,:)); % overidentification test stat 
                pvalET(r,:)=1-chi2cdf(J_statET(r,:),G-k); % pvalue of the test statistic
                if  pvalET(r,:)>=0.05 % at 5% sig level
                    c_ET(r,:)=1; % number of rejections
                end
                wald_ET(r,:)=wald_test(x,y,ng,G,betaET(r,:)',p_ET);
                pval_wald_ET(r,:)=1-chi2cdf(wald_ET(r,:),G-k);
                if  pval_wald_ET(r,:)>=0.05 % at 5% sig level
                    c_wald_ET(r,:)=1; % number of rejections
                end
                LM_ET(r,:)=LM_test(x,y,ng,G,betaET(r,:)',lambda,p_ET);
                pval_LM_ET(r,:)=1-chi2cdf(LM_ET(r,:),G-k);
                if  pval_LM_ET(r,:)>=0.05 % at 5% sig level
                    c_LM_ET(r,:)=1; % number of rejections
                end
            end %end of rep

            beta1='educ';
            b_names=char(strvcat(beta1));

            medbias2sls = median(beta2sls-onesr*beta_true');
            medbiasGMM = median(betaGMM-onesr*beta_true');

            medbiasEL = median(betaEL-onesr*beta_true');
            medbiasEL_biascorr = median(bias_crrctdEL-onesr*beta_true');

            medbiasET = median(betaET-onesr*beta_true');
            medbiasET_biascorr = median(bias_crrctdET-onesr*beta_true');

            covrateGMM=ones(1,k)-(onesr'*(abs((betaGMM-onesr*beta_true')./se_betaGMM)>1.6449*onesr*ones(1,k)))/rep;
            covrate2sls=ones(1,k)-(onesr'*(abs((beta2sls-onesr*beta_true')./se_beta2sls)>1.6449*onesr*ones(1,k)))/rep;

            covrateEL=ones(1,k)-(onesr'*(abs((betaEL-onesr*beta_true')./se_betaEL)>1.6449*onesr*ones(1,k)))/rep;
            covrateEL_biascorr=ones(1,k)-(onesr'*(abs((bias_crrctdEL-onesr*beta_true')./se_bias_crrctdEL)>1.6449*onesr*ones(1,k)))/rep; 

            covrateET=ones(1,k)-(onesr'*(abs((betaET-onesr*beta_true')./se_betaET)>1.6449*onesr*ones(1,k)))/rep;
            covrateET_biascorr=ones(1,k)-(onesr'*(abs((bias_crrctdET-onesr*beta_true')./se_bias_crrctdET)>1.6449*onesr*ones(1,k)))/rep; 

            fprintf('nrep    =  % 4.0f \n', rep );
            fprintf('Cohorts =  % 4.0f \n', C );
            fprintf('obs     =  % 4.0f \n', obs );
            fprintf('G       =  % 4.0f \n', G );
            fprintf('N       =  % 4.0f \n', mean(N));

            fprintf('k       =  % 4.0f \n', k );
            fprintf('beta    =  % s \n', b_names);
            fprintf('\n');
            fprintf('moments of u \n');
            fprintf('   mean        variance     third moment \n');
            for jj=1:G
                fprintf('% 4.6f    % 4.6f      % 4.6f    \n',mean(bar_ug(jj,:)),mean(sigma2_ug(jj,:)),mean(ug_3(jj,:)));
            end
            fprintf('\n');
            fprintf('Correlation of x and u \n');
            for jj=1:G
                for ii=1:k-1
                    fprintf('% 4.6f  \n',mean(xu_g(jj,ii,:)));
                end
            end
            fprintf('\n');
            fprintf('Estimators       Mean Estimate       Mean se     Median Estimate   Median se \n');
            for jj=1:k
            fprintf('2SLS              % 4.8f     % 4.8f     % 4.8f     % 4.8f\n', mean(beta2sls(:,jj)), mean(se_beta2sls(:,jj)), median(beta2sls(:,jj)), median(se_beta2sls(:,jj)));
            fprintf('GMM               % 4.8f     % 4.8f     % 4.8f     % 4.8f\n', mean(betaGMM(:,jj)), mean(se_betaGMM(:,jj)),median(betaGMM(:,jj)), median(se_betaGMM(:,jj)));
            fprintf('EL                % 4.8f     % 4.8f     % 4.8f     % 4.8f\n', mean(betaEL(:,jj)), mean(se_betaEL(:,jj)),median(betaEL(:,jj)), median(se_betaEL(:,jj)));
            fprintf('EL_biascorr       % 4.8f     % 4.8f     % 4.8f     % 4.8f\n', mean(bias_crrctdEL(:,jj)), mean(se_bias_crrctdEL(:,jj)), median(bias_crrctdEL(:,jj)), median(se_bias_crrctdEL(:,jj)));
            fprintf('ET                % 4.8f     % 4.8f     % 4.8f     % 4.8f\n', mean(betaET(:,jj)), mean(se_betaET(:,jj)),median(betaET(:,jj)), median(se_betaET(:,jj)));
            fprintf('ET_biascorr       % 4.8f     % 4.8f     % 4.8f     % 4.8f\n', mean(bias_crrctdET(:,jj)), mean(se_bias_crrctdET(:,jj)), median(bias_crrctdET(:,jj)), median(se_bias_crrctdET(:,jj)));
            end
            fprintf('\n');

            fprintf('Estimators         Median Bias    Coverage rate \n');
            for jj=1:k
            fprintf('2SLS               % 4.8f        % 4.4f\n', medbias2sls(:,jj), covrate2sls(:,jj));
            fprintf('GMM                % 4.8f        % 4.4f\n', medbiasGMM(:,jj), covrateGMM(:,jj));
            fprintf('EL                 % 4.8f        % 4.4f\n', medbiasEL(:,jj), covrateEL(:,jj));
            fprintf('EL_biascorr        % 4.8f        % 4.4f\n', medbiasEL_biascorr(:,jj), covrateEL_biascorr(:,jj));
            fprintf('ET                 % 4.8f        % 4.4f\n', medbiasET(:,jj), covrateET(:,jj));
            fprintf('ET_biascorr        % 4.8f        % 4.4f\n', medbiasET_biascorr(:,jj), covrateET_biascorr(:,jj));
            end
            fprintf('\n');
            fprintf('Number of acceptances of GMM  J test(5 percent sig level)        % 4.0f  \n', sum(c_GMM));
            fprintf('EL \n');

            fprintf('EL \n');
            fprintf('Number of acceptances of overid test(5 percent sig level)        % 4.0f  \n', sum(c_EL));
            fprintf('Number of acceptances of wald test(5 percent sig level)          % 4.0f  \n', sum(c_wald_EL));
            fprintf('Number of acceptances of LM test(5 percent sig level)            % 4.0f  \n', sum(c_LM_EL));
            fprintf('\n');
            fprintf('ET \n');
            fprintf('Number of acceptances of overid test(5 percent sig level)        % 4.0f  \n', sum(c_ET));
            fprintf('Number of acceptances of wald test(5 percent sig level)          % 4.0f  \n', sum(c_wald_ET));
            fprintf('Number of acceptances of LM test(5 percent sig level)            % 4.0f  \n', sum(c_LM_ET));

            fprintf('\n');
        end % end of different groups
    end % end of different obs
end % end of different cohorts

[Digitalsd]

In addition you can use the code:

stream = RandStream('mt19937ar','seed',1);
RandStream.setGlobalStream(stream);
randn

stream = RandStream('mt19937ar','seed',1);
RandStream.setGlobalStream(stream);
randn

"mt19937ar" indicates the used p.s.random generator algorithm