% calculations floaor the dog that doesn't bark. %revised version corresponds to paper revision, throws out lots of expermiments that got left out. See original "regressions.m" for those. clear all; close all; % good to start by clearing out stuff from previous runs % ************** % Switches and options % ***************** ntrial = 50000; % number of simulations for each one; Final uses 50000, but 5000 is PLENTY and gets to one decimal place print_graphs = 1; st = 1; % allows subsamples. 25 = 1950 doGW = 0; % do goyal welch? skip to save time if working on other stuff. redo_sim = 0; % redo simulation or load to play with graphs do_rx = 0; % do excess returns, "dividends" is the combination of d and rf do_real = 1; % choose ONE of real and rx. real is the base case; rx is a variation not pursued in the paper big_phi = 0; % do lots of phi? 0 will give data and 0.99 dobigrho = 0; % option to do rho one standard deviation higher. if ~((do_real==1)&(do_rx==0)&(st==1)&(dobigrho==0)); % this is the basic case print_graphs = 0; end; % phiv the vector of phis to try is selected below % ************* % Load data % ************** load all_vwr.txt; % load CRSP value weighted return with and without dividends cal = all_vwr(:,1); r = all_vwr(:,2); rx = all_vwr(:,3); load sbbi.txt; rf = sbbi(:,2); % 90 day bill cpi = sbbi(:,3); % inflation % -------------------------------------------- % 2. form dividend price ratio % D_t+1 (P_t+1 + D_t+1 ) P_t % --- --------------- ------ - 1 % P_t+1 P_t P_t+1 % % = R_t+1 / Rx_t+1 - 1 % then form dividend growth % D_t+1 / D_t = (D_t+1/P_t+1)/(D_t / P_t)*(P_t+1/P_t) % ----------------------------------------------------- dp = (1+r)./(1+rx) - 1; dd = [0 ; dp(2:end)./dp(1:end-1).*(1+rx(2:end))-1]; % do not use first dd cal = cal(st:end); dp = dp(st:end); dd = dd(st:end); r = r(st:end); rf = rf(st:end); cpi = cpi(st:end); rn = r; ddn = dd; rfn = rf; % preserves nominal versions if do_rx; disp(' *************** WARNING everything from here uses r = r/rf, dd = dd-rf'); r = (1+r)./(1+rf)-1; % use excess return in place of r -- excess LOG returns dd = (1+dd)./(1+rf)-1; % shifted rf from r to dd end; if do_real; disp(' *************** WARNING everything from here uses real returns and dd'); r = (1+r)./(1+cpi)-1; % use excess return in place of r rf = (1+rf)./(1+cpi)-1; % use excess return in place of r dd = (1+dd)./(1+cpi)-1; % shifted rf from r to dd end; % -------------------- % statstics % ------------------- disp('standard deviations'); fprintf(' %8.8s','sd(log(pd))','sd(logr)','sd(logdd)'); fprintf('\n'); fprintf(' %8.2f',std(log(dp)),std(log(1+r)),std(log(1+dd(2:end)))); fprintf('\n'); sddp_data = std(log(dp)); % save for later disp('covariance of pd, logr log dd') ; disp(cov([log(dp(2:end)) log(1+r(2:end)) log(1+dd(2:end))])); disp('correlation of pd, logr log dd') ; disp(corrcoef([log(dp(2:end)) log(1+r(2:end)) log(1+dd(2:end))])); disp('rho calculated from mean(dp), mean(pd) and exp(mean(log(p/d))') disp([ 1/mean(dp)/(1+1/mean(dp))... mean(1./dp)/(1+mean(1./dp))... exp(-mean(log(dp)))/(1+exp(-mean(log(dp))))]); PD = exp(-mean(log(dp))); rho = exp(-mean(log(dp)))/(1+exp(-mean(log(dp)))); disp('PD and rho'); disp([PD rho]); disp('mean log dp, sd dp, and standard error corrected for phi'); meandp = mean(log(dp)); disp(meandp); disp(std(log(dp))) T = size(dp,1); phi = [ones(T-1,1) log(dp(1:end-1))]\log(dp(2:end)); phi = phi(2); sedp = ( (var(log(dp))/(T-1)*(1+phi)/(1-phi))^0.5); disp('a higher rho -- mean log(dp) + sigma(log(dp) -- P/D and rho') PD2 = exp(-meandp+sedp); rho2 = PD2/(1+PD2); PD3 = exp(-meandp+2*sedp); rho3 = PD3/(1+PD3); disp('PD and rho +1 and +2 standard errors'); disp([PD2 rho2]); disp([PD3 rho3]); figure; plot((1926:2004)',-log(dp)) hold on plot((1926:2004)',-meandp*ones(T,1),'-r'); plot((1926:2004)',(-meandp+sedp)*ones(T,1),'--r'); plot((1926:2004)',(-meandp-sedp)*ones(T,1),'--r'); plot((1926:2004)',(-meandp+2*sedp)*ones(T,1),'--r'); plot((1926:2004)',(-meandp-2*sedp)*ones(T,1),'--r'); if dobigrho rho = rho2; PD = PD2; disp('********using higher PD and rho'); else; disp('Using rho from mean log dp'); end; fprintf('PD and rho %10.4f %10.4f \n', PD, rho); % ------------------------------------------- % regressions % ------------------------------------------- % 1. Table 1 regressions of returns on lagged returns disp('Forecasting regressions '); fprintf(' %10s %10s %10s %10s %10s \n','b','se','t(b)','R2','s(bx)'); disp(''); T = size(dp,1); % levels rhv = [ones(T-1,1) dp(1:T-1)]; lhv = r(2:end); [b,bse,r2] = olsgmm(lhv,rhv,0,0); bt = b./bse; fprintf('R on D/P %10.4f %10.4f %10.4f %10.4f %10.4f \n', [b(2) bse(2) bt(2) r2 std(rhv*b)]); rhv = [ones(T-1,1) dp(1:T-1)]; lhv = r(2:end)-rf(2:end); [b,bse,r2] = olsgmm(lhv,rhv,0,0); bt = b./bse; fprintf('R-Rf on D/P %10.4f %10.4f %10.4f %10.4f %10.4f \n', [b(2) bse(2) bt(2) r2 std(rhv*b)]); lhv = dd(2:end); [b,bse,r2] = olsgmm(lhv,rhv,0,0); bt = b./bse; fprintf('DD on D/P %10.4f %10.4f %10.4f %10.4f %10.4f\n', [b(2) bse(2) bt(2) r2 std(rhv*b)]); % logs rhv = [ones(T-1,1) log(dp(1:T-1))]; lhv = log(r(2:end)); [b,bse,r2] = olsgmm(lhv,rhv,0,0); bt = b./bse; fprintf('r on d-p %10.4f %10.4f %10.4f %10.4f %10.4f \n', [b(2) bse(2) bt(2) r2 std(rhv*b)]); lhv = log(dd(2:end)); [b,bse,r2] = olsgmm(lhv,rhv,0,0); bt = b./bse; fprintf('dd on d-p %10.4f %10.4f %10.4f %10.4f %10.4f\n', [b(2) bse(2) bt(2) r2 std(rhv*b)]); % not used in paper but hold in reserve rhv = [ones(T-5,1) dp(1:T-5)]; lhv = (1+r(2:end-4)).*(1+r(3:end-3)).*(1+r(4:end-2)).*(1+r(5:end-1)).*(1+r(6:end))-1; [b,bse,r2] = olsgmm(lhv,rhv,0,0); bt = b./bse; fprintf('R5 on D/P %10.4f %10.4f %10.4f %10.4f %10.4f \n', [b(2) bse(2) bt(2) r2 std(rhv*b)]); rhv = [ones(T-10,1) dp(1:T-10)]; lhv = (1+r(2:end-9)).*(1+r(3:end-8)).*(1+r(4:end-7)).*(1+r(5:end-6)).*(1+r(6:end-5)).*(1+r(7:end-4))... .*(1+r(8:end-3)).*(1+r(9:end-2)).*(1+r(10:end-1)).*(1+r(11:end))-1; [b,bse,r2] = olsgmm(lhv,rhv,0,0); bt = b./bse; fprintf('R10 on D/P %10.4f %10.4f %10.4f %10.4f %10.4f \n', [b(2) bse(2) bt(2) r2 std(rhv*b)]); % ********************************* % Main two-variable VAR regressions % ********************************* rhv = [ones(T-1,1) log(dp(1:T-1))]; lhv = log(1+r(2:end)); [b,bse,r2] = olsgmm(lhv,rhv,0,0); bt = b./bse; brdp_data = b(2); trdp_data = bt(2); %save for comparision on monte carlos err_r = lhv - rhv*b; fprintf('log r on log dp %10.4f %10.4f %10.4f %10.4f %10.4f\n', [b(2) bse(2) bt(2) r2 std(rhv*b)]); rhv = [ones(T-1,1) log(dp(1:T-1))]; lhv = log(1+dd(2:end)); [b,bse,r2] = olsgmm(lhv,rhv,0,0); bt = b./bse; bddp_data = b(2); tddp_data = bt(2); %save for comparision on monte carlos err_d = lhv - rhv*b; fprintf('log dd on log dp %10.4f %10.4f %10.4f %10.4f %10.4f\n', [b(2) bse(2) bt(2) r2 std(rhv*b)]); rhv = [ones(T-1,1) log(dp(1:T-1))]; lhv = log(dp(2:end)); [b,bse,r2] = olsgmm(lhv,rhv,0,0); bt = b./bse; phidp_data = b(2); tphidp_data = bt(2); %save for comparision on monte carlos err_dp = lhv - rhv*b; fprintf('log dp on log dp %10.4f %10.4f %10.4f %10.4f %10.4f\n', [b(2) bse(2) bt(2) r2 std(rhv*b)]); disp('Error sd/correlation matrix'); fprintf(' %10s %10s %10s \n','r','d','dp'); fprintf(' r %10.3f %10.3f %10.3f \n',std(err_r),corr(err_r,err_d),corr(err_r,err_dp)); fprintf(' dd %10.3f %10.3f %10.3f \n',corr(err_r,err_d),std(err_d),corr(err_d,err_dp)); fprintf(' dp %10.3f %10.3f %10.3f \n',corr(err_r,err_dp),corr(err_d,err_dp),std(err_dp)); var_err_dp_data = var(err_dp); % used later % Implied estimates : identity is not exact, so use two regressions, imply the third coefficient brdp_data_impl = 1 - rho*phidp_data + bddp_data; bddp_data_impl = rho*phidp_data - 1 + brdp_data; phidp_data_impl = 1/rho*(bddp_data - brdp_data +1); disp('implied vs direct coefficients'); fprintf(' %10s %10s %10s \n','b_r','b_d','phi'); fprintf(' direct %10.3f %10.3f %10.3f \n',brdp_data,bddp_data,phidp_data); fprintf(' implied %10.3f %10.3f %10.3f \n',brdp_data_impl,bddp_data_impl,phidp_data_impl); err_r_impl = err_d - rho*err_dp; err_d_impl = err_r + rho*err_dp; err_dp_impl = (err_d - err_r)/rho; disp('implied error sd/correlation matrix -- implied error on rows'); fprintf(' %10s %10s %10s \n','r','d','dp'); fprintf(' r %10.3f %10.3f %10.3f \n',std(err_r_impl),corr(err_r_impl,err_d),corr(err_r_impl,err_dp)); fprintf(' dd %10.3f %10.3f %10.3f \n',corr(err_r,err_d_impl),std(err_d_impl),corr(err_d_impl,err_dp)); fprintf(' dp %10.3f %10.3f %10.3f \n',corr(err_r,err_dp_impl),corr(err_d,err_dp_impl),std(err_dp_impl)); % ****************** % "Long run" regressions, i.e. estimate b_r/(1-rho phi) by estimating br and phi % in each case use one two-variable VAR so the identity works % ****************** % long run estimates using return - dd var brlr = brdp_data/(1-rho*phidp_data_impl); bdlr = bddp_data/(1-rho*phidp_data_impl); cov_r_dp = cov(err_r,err_dp_impl); cov_r_dp = cov_r_dp(1,2); cov_d_dp = cov(err_d,err_dp_impl); cov_d_dp = cov_d_dp(1,2); sebrlr = 1/(1-rho*phidp_data_impl)*(var(err_r) + 2*rho*brlr*cov_r_dp + rho^2*brlr^2*var(err_dp_impl))^0.5/((T-1)^0.5*std(log(dp(1:T-1)))); sebdlr = 1/(1-rho*phidp_data_impl)*(var(err_d) + 2*rho*bdlr*cov_d_dp + rho^2*bdlr^2*var(err_dp_impl))^0.5/((T-1)^0.5*std(log(dp(1:T-1)))); disp(''); disp('long-run regressions -- r, dd VAR, dp implied') fprintf(' %10s %10s %10s \n','b','se','t'); fprintf(' r %10.3f %10.3f %10.3f \n', brlr, sebrlr, brlr/sebrlr); fprintf(' d %10.3f %10.3f %10.3f \n', bdlr, sebdlr, (bdlr+1)/sebdlr); % long run estimates using dd - dp var brlr = brdp_data_impl/(1-rho*phidp_data); bdlr = bddp_data/(1-rho*phidp_data); cov_r_dp = cov(err_r_impl,err_dp); cov_r_dp = cov_r_dp(1,2); cov_d_dp = cov(err_d,err_dp); cov_d_dp = cov_d_dp(1,2); sebrlr = 1/(1-rho*phidp_data)*(var(err_r_impl) + 2*rho*brlr*cov_r_dp + rho^2*brlr^2*var(err_dp))^0.5/((T-1)^0.5*std(log(dp(1:T-1)))); sebdlr = 1/(1-rho*phidp_data)*(var(err_d) + 2*rho*bdlr*cov_d_dp + rho^2*bdlr^2*var(err_dp))^0.5/((T-1)^0.5*std(log(dp(1:T-1)))); disp(''); disp('long-run regressions -- dd - dp VAR, r implied') fprintf(' %10s %10s %10s \n','b','se','t'); fprintf(' r %10.3f %10.3f %10.3f \n', brlr, sebrlr, brlr/sebrlr); fprintf(' d %10.3f %10.3f %10.3f \n', bdlr, sebdlr, (bdlr+1)/sebdlr); % choose one for simulation! brlr_data = brlr; bdlr_data = bdlr; trlr_data = brlr/sebrlr; tdlr_data = (bdlr-(-1))/sebdlr; % long run estimates using return - dp , dd implied brlr = brdp_data/(1-rho*phidp_data); bdlr = bddp_data_impl/(1-rho*phidp_data); cov_r_dp = cov(err_r,err_dp); cov_r_dp = cov_r_dp(1,2); cov_d_dp = cov(err_d_impl,err_dp); cov_d_dp = cov_d_dp(1,2); sebrlr = 1/(1-rho*phidp_data)*(var(err_r) + 2*rho*brlr*cov_r_dp + rho^2*brlr^2*var(err_dp))^0.5/((T-1)^0.5*std(log(dp(1:T-1)))); sebdlr = 1/(1-rho*phidp_data)*(var(err_d_impl) + 2*rho*bdlr*cov_d_dp + rho^2*bdlr^2*var(err_dp))^0.5/((T-1)^0.5*std(log(dp(1:T-1)))); disp(''); disp('long-run regressions -- r, dp VAR, dd implied') fprintf(' %10s %10s %10s \n','b','se','t'); fprintf(' r %10.3f %10.3f %10.3f \n', brlr, sebrlr, brlr/sebrlr); fprintf(' d %10.3f %10.3f %10.3f \n', bdlr, sebdlr, (bdlr+1)/sebdlr); disp(''); % **************** % univariate likelihood of phi % ******************** phix = (0.7:0.001:0.999)'; phix1 = 1+0*phix; sumt = log(dp(2:T))*phix1'-log(dp(1:T-1))*phix'; sumt = sum(sumt)'; % now a vertical vector conformable with phix with the sum term for each phi a = 1./(T + 2*phix./(1-phix)) .* ( log(dp(1))*(1+phix) + sumt); sumt = (log(dp(2:T))*phix1'-ones(T-1,1)*a'-log(dp(1:T-1))*phix').^2; sumt = sum(sumt)'; % now a vertical vector conformable with phix with the sum term for each phi sig2 = 1/T*( (log(dp(1))-a./(1-phix)).^2.*(1-phix).^2 + sumt); Lu = -T/2*log(2*pi) - 1/2*log(sig2./(1-phix.^2))... - 1/2.*((log(dp(1))-a./(1-phix)).^2).*(1-phix.^2)./sig2... -(T-1)/2*log(sig2) - 1/2*sumt./sig2; [Lumax, indLu] = max(Lu); % conditional phixc = (0.7:0.001:1.05)'; phixc1 = 1+0*phixc; sumt = log(dp(2:T))*phixc1'-log(dp(1:T-1))*phixc'; sumt = sum(sumt)'; % now a vertical vector conformable with phix with the sum term for each phi a = 1./(T-1).*sumt; sumt = (log(dp(2:T))*phixc1'-ones(T-1,1)*a'-log(dp(1:T-1))*phixc').^2; sumt = sum(sumt)'; % now a vertical vector conformable with phix with the sum term for each phi sig2 = 1/(T-1) * sumt; Lc = -(T-1)/2*log(2*pi) -(T-1)/2*log(sig2) - 1/2*sumt./sig2; [Lcmax, indLc] = max(Lc); fprintf(' MaxL of AR(1) log dp %10s %10s %10s %10s \n', 'Lu', 'phiu', 'Lc','phic'); fprintf(' %10.3f %10.3f %10.3f %10.3f \n',Lumax,phix(indLu),Lcmax,phix(indLc)); figure; plot(phix,Lu,'-b'); hold on; plot(phixc,Lc,'-g') plot(phix(indLu),Lumax,'vb'); plot(phixc(indLc),Lcmax,'og'); xlabel('\phi'); ylabel('Log L'); legend('Unconditional','Conditional'); figure; plot([phix;1],[exp(Lu)/sum(exp(Lu));0],'-b','LineWidth',2); hold on; plot(phixc,exp(Lc)/sum(exp(Lc)),'--r','LineWidth',2) plot(phix(indLu),exp(Lumax)/sum(exp(Lu)),'vb'); plot(phixc(indLc),exp(Lcmax)/sum(exp(Lc)),'or'); xlabel('\phi'); ylabel('Likelihood'); axis([0.8 1.05 0 0.012]); text(phix(indLu),exp(Lumax)/sum(exp(Lu))+0.0005,'Uconditional','Color','b','BackgroundColor',[1 1 1]) %text(phixc(indLc),exp(Lcmax)/sum(exp(Lc))+0.0005,'Conditional','Color','r','BackgroundColor',[1 1 1]) indx = find(phixc==1); text(phixc(indx),exp(Lc(indx))/sum(exp(Lc))+0.0005,'Conditional','Color','r','BackgroundColor',[1 1 1]) %plot(phix(indLu)*[1 1]',[0 exp(Lumax)]','--k') %plot(phixc(indLc)*[1 1]',[0 exp(Lcmax)]','--k') %legend('Unconditional','Conditional'); set(gca,'Ytick',[]); figheight =4; figwidth = 6; % set here common size for all figures in inches. This is for two panel pictures set(gcf, 'PaperPositionMode', 'manual'); set(gcf, 'PaperUnits', 'inches'); set(gcf, 'PaperPosition', [0.25 2.5 figwidth figheight]); if print_graphs; print -depsc2 phi_likelihood.eps; end; % cumulative for sampling FLu = cumsum(exp(Lu)); FLu = FLu/FLu(end); % ************************** % long-horizon regressions etc. in logs for long-horizon section % ************************** % delta d t+j on dpt hmax = 25; bvd = zeros(hmax,1); bsevd = zeros(hmax,1); bvr = zeros(hmax,1); bsevr = zeros(hmax,1); for h = 1:hmax; rhv = [ones(T-h,1) log(dp(1:T-h))]; lhv = log(1+dd(1+h:end)); [b,bse,r2] = olsgmm(lhv,rhv,2*h,1); bvd(h) = b(2); bsevd(h) = bse(2); lhv = log(1+r(1+h:end)); [b,bse,r2] = olsgmm(lhv,rhv,2*h,1); [b2,bse2,r2] = olsgmm(lhv,rhv,h,2); bvr(h) = b(2); bsevr(h) = bse(2); end; figure; subplot(2,1,1); plot((1:hmax)',bvd,'-v','LineWidth',2) hold on plot((1:hmax)',bsevd*[-2 2],':k','LineWidth',2) plot((1:hmax)',0*bvd,'-k') plot((1:hmax)',bddp_data*phidp_data.^(0:hmax-1)','--r','LineWidth',2); axis([0 hmax -0.15 0.15]); title('Dividend growth'); subplot(2,1,2); plot((1:hmax)',bvr,'-v','LineWidth',2) hold on plot((1:hmax)',bsevr*[-2 2],':k','LineWidth',2) plot((1:hmax)',0*bvr,'-k') plot((1:hmax)',brdp_data*phidp_data.^(0:hmax-1)','--r','LineWidth',2); axis([0 hmax -0.15 0.15]); title('Returns'); xlabel('Horizon in years'); if print_graphs; print -depsc2 long_horizon_1.eps; end; % using sums of the single regressions rather htan the direct regression of % multiyear returns impl_lr_d = cumsum(rho.^(0:hmax-1)'.*bvd); impl_lr_r = cumsum(rho.^(0:hmax-1)'.*bvr); figure; subplot(2,1,1); plot((1:hmax)',cumsum(rho.^(0:hmax-1)'.*bvd),'-v','LineWidth',2) hold on plot((1:hmax)',0*bvd,'-k') plot((1:hmax)',bddp_data*(1-(rho*phidp_data).^(1:hmax)')./(1-rho*phidp_data),'--r','LineWidth',2); axis([0 hmax -0.3 1.2]); title('Dividend growth'); subplot(2,1,2); plot((1:hmax)',cumsum(rho.^(0:hmax-1)'.*bvr),'-v','LineWidth',2) hold on plot((1:hmax)',0*bvr,'-k') plot((1:hmax)',brdp_data*(1-rho*phidp_data.^(1:hmax)')./(1-rho*phidp_data),'--r','LineWidth',2); axis([0 hmax -0.3 1.2]); title('Returns'); xlabel('Horizon in years'); %print -depsc2 long_horizon_2.eps % direct long horizon estimates % sum rho^ j d_t+j on dp t % hmax = 25; bvd = zeros(hmax,1); bsevd = zeros(hmax,1); bsevd2 = zeros(hmax,1); bvr = zeros(hmax,1); bvr2 = zeros(hmax,1); bsevr = zeros(hmax,1); bsevr2 = zeros(hmax,1); for h = 1:hmax; rhv = [ones(T-h,1) log(dp(1:T-h))]; lhv = log(1+dd(1+1:T-h+1)); for j = 2:h lhv = lhv+rho^(j-1)*log(1+dd(1+j:T-h+j)); end; [b,bse,r2] = olsgmm(lhv,rhv,2*h,1); [b2,bse2,r2] = olsgmm(lhv,rhv,h,rho); bvd(h) = b(2); bsevd(h) = bse(2); bsevd2(h) = bse2(2); lhv = log(1+r(1+1:T-h+1)); for j = 2:h lhv = lhv+rho^(j-1)*log(1+r(1+j:T-h+j)); end; [b,bse,r2] = olsgmm(lhv,rhv,2*h,1); [b2,bse2,r2] = olsgmm(lhv,rhv,h,rho); bvr(h) = b(2); bvr2(h) = b2(2); bsevr(h) = bse(2); bsevr2(h) = bse2(2); end; lr1_data = bvr(1); lr5_data = bvr(5); lr10_data = bvr(10); lr15_data = bvr(15); % save to compare with monte carlo lr20_data = bvr(20); ld1_data = bvd(1); ld5_data = bvd(5); ld10_data = bvd(10); ld15_data = bvd(15); % save to compare with monte carlo ld20_data = bvd(20); figure; subplot(2,1,1); plot((1:hmax)',bvd,'-v','LineWidth',2) hold on plot((1:hmax)',impl_lr_d,'-ro','LineWidth',2) %plot((1:hmax)',bsevd*[-2 2],':k','LineWidth',2) %plot((1:hmax)',bsevd2*[-2 2],':k','LineWidth',2) plot((1:hmax)',0*bvd,'-k') plot((1:hmax)',bddp_data*(1-rho*phidp_data.^(1:hmax)')./(1-rho*phidp_data),'--m','LineWidth',2); axis([0 hmax -0.2 1.55]); title('Dividend growth'); subplot(2,1,2); plot((1:hmax)',bvr,'-v','LineWidth',2) hold on plot((1:hmax)',impl_lr_r,'-ro','LineWidth',2) %plot((1:hmax)',bsevr*[-2 2],':k','LineWidth',2) %plot((1:hmax)',bsevr2*[-2 2],':k','LineWidth',2) plot((1:hmax)',0*bvr,'-k') plot((1:hmax)',brdp_data*(1-rho*phidp_data.^(1:hmax)')./(1-rho*phidp_data),'--m','LineWidth',2); axis([0 hmax -0.2 1.55]); title('Returns'); xlabel('Horizon in years'); if print_graphs print -depsc2 long_horizon_2.eps end; % direct long horizon estimates % unweighted long horizon returns sum d_t+j on dp t % hmax = 25; bvdu = zeros(hmax,1); bsevdu = zeros(hmax,1); bsevd2u = zeros(hmax,1); bvru = zeros(hmax,1); bvr2u = zeros(hmax,1); bsevru = zeros(hmax,1); bsevr2u = zeros(hmax,1); for h = 1:hmax; rhv = [ones(T-h,1) log(dp(1:T-h))]; lhv = log(1+dd(1+1:T-h+1)); for j = 2:h lhv = lhv+log(1+dd(1+j:T-h+j)); end; [b,bse,r2] = olsgmm(lhv,rhv,2*h,1); [b2,bse2,r2] = olsgmm(lhv,rhv,h,rho); bvdu(h) = b(2); bsevdu(h) = bse(2); bsevd2u(h) = bse2(2); lhv = log(1+r(1+1:T-h+1)); for j = 2:h lhv = lhv+log(1+r(1+j:T-h+j)); end; [b,bse,r2] = olsgmm(lhv,rhv,2*h,1); [b2,bse2,r2] = olsgmm(lhv,rhv,h,rho); bvru(h) = b(2); bvr2u(h) = b2(2); bsevru(h) = bse(2); bsevr2u(h) = bse2(2); end; lru1_data = bvru(1); lru5_data = bvru(5); lru10_data = bvru(10); lru15_data = bvru(15); % save to compare with monte carlo lru20_data = bvru(20); ldu1_data = bvdu(1); ldu5_data = bvdu(5); ldu10_data = bvdu(10); ldu15_data = bvdu(15); % save to compare with monte carlo ldu20_data = bvdu(20); % ********************************************************************** % find Goyal-welch statistic -- also used to compare simulations below % ********************************************************************* if st == 1; % not reprogrammed for shorter sample -- silly to do it. err = zeros(T,1); err2 = zeros(T,1); % starting at t=20, forecast returns using data up to t for i = 20:T-1; [bsim,bsesim,r2sim] = olsgmm(log(1+r(2:i)),[ones(i-1,1) log(dp(1:i-1))],0,0); fcstr = [1 log(dp(i))]*bsim; err(i+1)=log(1+r(i+1))-fcstr; % forecast using dp fcst2 = mean(log(1+r(1:i))); err2(i+1) = log(1+r(i+1))-fcst2; % forecast using mean end; err = err(21:end); err2 = err2(21:end); figure; plot((1946:2004)',-[(cumsum(err.^2).^0.5)-(cumsum(err2.^2).^0.5)]); title('rmse of dp - rmse of constant'); disp('rmse of dp, r2 of mean, difference (const-dp)'); msedp = mean(err.^2)^0.5; msecon = mean(err2.^2)^0.5 ; disp([msedp msecon msecon-msedp]); % Goyal-Welch at longer horizons h = 5; % return-forecast horizon in years err = zeros(T,1); err2 = zeros(T,1); fcstrt = zeros(T,1); fcstrcont = zeros(T,1); bt = zeros(T,1); % starting at t=20, forecast returns using data up to t lhr = filter(ones(h,1),1,log(1+r)); % full sample forecast [blh,bselh,r2lh] = olsgmm(lhr(1+h:T),[ones(T-h,1) log(dp(1:T-h))],0,0); fcstr_full = [ones(T,1) log(dp(1:T))]*blh; for i = 20:T-h; [bsimlh,bsesimlh,r2simlh] = olsgmm(lhr(1+h:i),[ones(i-h,1) log(dp(1:i-h))],0,0); fcstr = [1 log(dp(i))]*bsimlh; bt(i) = bsimlh(2); fcstrt(i) = fcstr; err(i+h)=lhr(i+h)-fcstr; % forecast using dp fcst2 = mean(lhr(h:i)); fcstrcont(i) = fcst2; err2(i+h) = lhr(i+h)-fcst2; % forecast using mean end; err = err(20+h:end); err2 = err2(20+h:end); figure; plot((1926:2004-h)',[fcstr_full(1:T-h) fcstrt(1:T-h) fcstrcont(1:T-h) lhr(1+h:T)]); xlabel('Date forecast made'); legend('full sample','ex ante','constant','actual'); title(['forecast and actual ' int2str(h) ' year returns']); figure; plot((1926:2004-h)',[bt(1:T-h) ones(T-h,1)*blh(2)]); xlabel('Date forecast made'); legend('ex ante','full sample'); title('regression coefficient on dp'); figure; plot((1945+h:2004)',-[(cumsum(err.^2)).^0.5-(cumsum(err2.^2)).^0.5]); title('Long horizon rmse of constant - rmse of dp'); disp('rmse of dp, r2 of mean, difference (const-dp)'); msedplh = mean(err.^2)^0.5; mseconlh = mean(err2.^2)^0.5; disp([msedplh mseconlh mseconlh-msedplh]); end; % ends if st == 1; % ************************** % simulations for Er=const % *************************** phiv = [phidp_data 0.99]'; % -99 codes for picking phi if big_phi phiv = [0.90 phidp_data 0.96 0.98 0.99 1.0 1.01 -99 ]'; end; lotsprint = 1; % print big table or not below if size(phiv) > 2; lotsprint = 0 ; end; N = size(phiv,1); brdp_sim = zeros(ntrial,N); bddp_sim= zeros(ntrial,N); bdpdp_sim= zeros(ntrial,N); trdp_sim = zeros(ntrial,N); tddp_sim = zeros(ntrial,N); tdpdp_sim = zeros(ntrial,N); sddp_sim = zeros(ntrial,N); phi_sim = zeros(ntrial,N); brlr_sim = zeros(ntrial,N); bdlr_sim = zeros(ntrial,N); trlr_sim = zeros(ntrial,N); tdlr_sim = zeros(ntrial,N); phinull_sim = zeros(ntrial,N); lr1 = zeros(ntrial,N); lr5 = zeros(ntrial,N); lr10 = zeros(ntrial,N); lr15 = zeros(ntrial,N); lr20 = zeros(ntrial,N); ld1 = zeros(ntrial,N); ld5 = zeros(ntrial,N); ld10 = zeros(ntrial,N); ld15 = zeros(ntrial,N); ld20 = zeros(ntrial,N); lru1 = zeros(ntrial,N); lru5 = zeros(ntrial,N); lru10 = zeros(ntrial,N); lru15 = zeros(ntrial,N); lru20 = zeros(ntrial,N); ldu1 = zeros(ntrial,N); ldu5 = zeros(ntrial,N); ldu10 = zeros(ntrial,N); ldu15 = zeros(ntrial,N); ldu20 = zeros(ntrial,N); msediff = zeros(ntrial,N); randn('state',0); % intitialize so graphs look the same covshock = cov(err_dp,err_d); % Simulatin uses dp - d var and implies the return sigma_dp = covshock(1,1)^0.5; sigma_d = covshock(2,2)^0.5; sigma_d_dp = covshock(1,2); if redo_sim; for phitry = 1:size(phiv,1); disp('************'); phi = phiv(phitry); % 0.94; % fix phi if phi == -99; draw_phi = 1; disp('drawing phi each time from likelihood'); else; draw_phi = 0; fprintf(' assumed phi value and corresponding bd values \n'); fprintf(' %10.4f %10.4f \n',phi, 1-rho*phi); end; for trial = 1:ntrial; if draw_phi == 1; % codes to pick phi from likelihood p = rand; [prob,indx] = min(abs(FLu-p)); phi = phix(indx); end; phinull_sim(trial,phitry) = phi; % keep track of what phi you're using e_dp_t = sigma_dp*randn(T,1); % y = cov(xy)/var(x) x + e, % var(y) = var(x)*cov(xy)^2/var(x)^2 + var(e) % var(y) -cov(xy)^2/var(x) = var(e) e_d_t = sigma_d_dp/sigma_dp^2*e_dp_t + (sigma_d^2 - sigma_d_dp^2/sigma_dp^2)^0.5*randn(T,1); % pick initial dp from its unconditional density if phi < 1; dp0 = sigma_dp/(1-phi^2)^0.5*randn; else; dp0 = 0; end; dpt = filter(1,[1 -phi],[dp0; e_dp_t]); dpt = dpt(2:end); ddt = [NaN; (rho*phi-1)*dpt(1:T-1)] + e_d_t; % first dd is never used, so no worries here rt = [ NaN; ddt(2:T) + dpt(1:T-1) - rho*dpt(2:T)]; % regression of r, dd on dp rhv = [ones(T-1,1) dpt(1:T-1) ]; lhv = [rt(2:end) ddt(2:end) dpt(2:end)]; [b,bse,r2] = olsgmm(lhv,rhv,0,0); errs = lhv-rhv*b; err_r = errs(:,1); err_dd = errs(:,2); err_dp = errs(:,3); bt = b./bse; % keep track of things across trials; brdp_sim(trial,phitry) = b(2,1); bddp_sim(trial,phitry) = b(2,2); bdpdp_sim(trial,phitry) = b(2,3); brdp = b(2,1); bddp = b(2,2); bdpdp = b(2,3); trdp_sim(trial,phitry) = bt(2,1); tddp_sim(trial,phitry) = bt(2,2); tdpdp_sim(trial,phitry) = bt(2,3); sddp_sim(trial,phitry) = std(dpt); % implied "long run" estimates brlr = brdp/(1-rho*bdpdp); bdlr = bddp/(1-rho*bdpdp); brlr_sim(trial,phitry) = brlr; bdlr_sim(trial,phitry) = bdlr; cov_r_dp = cov(err_r,err_dp); cov_r_dp = cov_r_dp(1,2); cov_d_dp = cov(err_d,err_dp); cov_d_dp = cov_d_dp(1,2); sebrlr = 1/(1-rho*bdpdp)*(var(err_r) + 2*rho*brlr*cov_r_dp + rho^2*brlr^2*var(err_dp))^0.5/... ((T-1)^0.5*std(dpt(1:T-1))); sebdlr = 1/(1-rho*bdpdp)*(var(err_d) + 2*rho*bdlr*cov_d_dp + rho^2*bdlr^2*var(err_dp))^0.5/... ((T-1)^0.5*std(dpt(1:T-1))); trlr_sim(trial,phitry) = brlr/sebrlr; tdlr_sim(trial,phitry) = (bdlr-(-1))/sebdlr; % direct long run estimates for h = [ 1 5 10 15 20]; rhv = [ones(T-h,1) dpt(1:T-h)]; lhv = ddt(1+1:T-h+1); for j = 2:h lhv = lhv+rho^(j-1)*ddt(1+j:T-h+j); end; [bd,bse,r2] = olsgmm(lhv,rhv,2*h,1); %[b2,bse2,r2] = olsgmm(lhv,rhv,h,rho); lhv = rt(1+1:T-h+1); for j = 2:h lhv = lhv+rho^(j-1)*rt(1+j:T-h+j); end; [b,bse,r2] = olsgmm(lhv,rhv,2*h,1); %[b2,bse2,r2] = olsgmm(lhv,rhv,h,rho); if h == 1; lr1(trial,phitry) = b(2); ld1(trial,phitry) = bd(2); elseif h == 5; lr5(trial,phitry) = b(2); ld5(trial,phitry) = bd(2); elseif h == 10; lr10(trial,phitry) = b(2); ld10(trial,phitry) = bd(2); elseif h == 15; lr15(trial,phitry) = b(2); ld15(trial,phitry) = bd(2); elseif h == 20; lr20(trial,phitry) = b(2); ld20(trial,phitry) = bd(2); end; end; % same unweighted for h = [ 1 5 10 15 20]; rhv = [ones(T-h,1) dpt(1:T-h)]; lhv = ddt(1+1:T-h+1); for j = 2:h lhv = lhv+ddt(1+j:T-h+j); end; [bd,bse,r2] = olsgmm(lhv,rhv,2*h,1); %[b2,bse2,r2] = olsgmm(lhv,rhv,h,rho); lhv = rt(1+1:T-h+1); for j = 2:h lhv = lhv+rt(1+j:T-h+j); end; [b,bse,r2] = olsgmm(lhv,rhv,2*h,1); %[b2,bse2,r2] = olsgmm(lhv,rhv,h,rho); if h == 1; lru1(trial,phitry) = b(2); ldu1(trial,phitry) = bd(2); elseif h == 5; lru5(trial,phitry) = b(2); ldu5(trial,phitry) = bd(2); elseif h == 10; lru10(trial,phitry) = b(2); ldu10(trial,phitry) = bd(2); elseif h == 15; lru15(trial,phitry) = b(2); ldu15(trial,phitry) = bd(2); elseif h == 20; lru20(trial,phitry) = b(2); ldu20(trial,phitry) = bd(2); end; end; % Goyal welch if doGW; % allows you to skip and save time err = zeros(T,1); err2 = zeros(T,1); % starting at t=20, forecast returns using data up to t for i = 20:T-1; [bsim,bsesim,r2sim] = olsgmm(rt(2:i),[ones(i-1,1) dpt(1:i-1)],0,0); fcstr = [1 dpt(i)]*bsim; err(i+1)= rt(i+1)-fcstr; % forecast using dp fcst2 = mean(rt(2:i)); err2(i+1) = rt(i+1)-fcst2; % forecast using mean end; err = err(21:end); err2 = err2(21:end); msedpsim = mean(err.^2)^0.5; mseconsim = mean(err2.^2)^0.5; msediff(trial,phitry) = mseconsim-msedpsim; end; end; % ends trials end; %ends phi end; % ends if redo_sim if dobigrho == 1; % NOTE this option does not check other options are set right -- be careful if redo_sim == 1; save simresults_dobigrho else load simresults_dobigrho end end if dobigrho == 0; if (ntrial == 50000) & (redo_sim == 1) & (do_real == 1) & (big_phi == 0 ) save simresults_real_twophi elseif (ntrial == 50000) & (redo_sim == 0) & (do_real == 1) & (big_phi == 0 ) load simresults_real_twophi redo_sim = 0; else; if doGW & redo_sim; save simresults; % save the whole thing so you cna work on figures without rerunning. end; if ~doGW & redo_sim; save simresults_noGW; end; if doGW & (~redo_sim); load simresults; % save the whole thing so you cna work on figures without rerunning. end; if ~doGW & (~redo_sim); load simresults_noGW; end; end; end; % big table of results if lotsprint; for phitry = 1:size(phiv,1); disp(['******* results for phi = ' num2str(phiv(phitry),'%5.2f')]); % give statistics -- graphs below fprintf(' probability of %10s %10s %10s %10s %10s \n','br > data','bdd > data','sdpbrdp_data)/ntrial,... 100*sum(bddp_sim(:,phitry)>bddp_data)/ntrial,... 100*sum(sddp_sim(:,phitry)brdp_data)&(bddp_sim(:,phitry)>bddp_data))/ntrial,... 100*sum((brdp_sim(:,phitry)>brdp_data)&(bddp_sim(:,phitry)>bddp_data)&(sddp_sim(:,phitry) data','phi&bdd','phi&br'); fprintf(' %10.5f %10.5f %10.5f \n',... 100*sum(bdpdp_sim(:,phitry)>phidp_data)/ntrial,... 100*sum((bdpdp_sim(:,phitry)>phidp_data)&(bddp_sim(:,phitry)>bddp_data))/ntrial,... 100*sum((bdpdp_sim(:,phitry)>phidp_data)&(brdp_sim(:,phitry)>brdp_data))/ntrial); fprintf(' same, phi id %10s %10s %10s \n','phi > data','phi&bdd','phi&br'); fprintf(' %10.5f %10.5f %10.5f \n',... 100*sum(bdpdp_sim(:,phitry)>phidp_data_impl)/ntrial,... 100*sum((bdpdp_sim(:,phitry)>phidp_data_impl)&(bddp_sim(:,phitry)>bddp_data))/ntrial,... 100*sum((bdpdp_sim(:,phitry)>phidp_data_impl)&(brdp_sim(:,phitry)>brdp_data))/ntrial); fprintf(' probability of %10s %10s \n','phi&bd&br','phi_id&&bd&br'); fprintf(' %10.5f %10.5f \n',... 100*sum((bdpdp_sim(:,phitry)>phidp_data)&(bddp_sim(:,phitry)>bddp_data)&(brdp_sim(:,phitry)>brdp_data))/ntrial,... 100*sum((bdpdp_sim(:,phitry)>phidp_data_impl)&(bddp_sim(:,phitry)>bddp_data)&(brdp_sim(:,phitry)>brdp_data))/ntrial); fprintf(' probability of %10s %10s %10s \n','tphi > data','tphi&tdd','tphi&tr'); fprintf(' %10.5f %10.5f %10.5f \n',... 100*sum(tdpdp_sim(:,phitry)>tphidp_data)/ntrial,... 100*sum((tdpdp_sim(:,phitry)>tphidp_data)&(tddp_sim(:,phitry)>tddp_data))/ntrial,... 100*sum((tdpdp_sim(:,phitry)>tphidp_data)&(trdp_sim(:,phitry)>trdp_data))/ntrial); fprintf(' probability of %10s %10s %10s \n','tr > tdata','tdd > tdata','tboth'); fprintf(' %10.5f %10.5f %10.5f \n', ... 100*sum(trdp_sim(:,phitry)>trdp_data)/ntrial,... 100*sum(tddp_sim(:,phitry)>tddp_data)/ntrial,... 100*sum((trdp_sim(:,phitry)>trdp_data)&(tddp_sim(:,phitry)>tddp_data))/ntrial); fprintf(' means %10s %10s %10s %10s %10s %10s \n','br','bd','phi','brlr','bdlr'); fprintf(' %10.3f %10.3f %10.3f %10.3f %10.3f \n',... mean(brdp_sim(:,phitry)),mean(bddp_sim(:,phitry)),mean(bdpdp_sim(:,phitry)),mean(brlr_sim(:,phitry)),mean(bdlr_sim(:,phitry))); disp('Statistics for long run forecasts'); fprintf(' %10s %10s %10s %10s %10s %10s %10s %10s %10s\n', 'phi', 'br','bd','br,impl br','bd,impl bd','br,impl phi','bd,impl phi','E(br)','E(bd)'); fprintf(' %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f\n',... phiv(phitry),... 100*sum(brdp_sim(:,phitry)./(1-rho*bdpdp_sim(:,phitry)) > brdp_data/(1-rho*phidp_data))/ntrial,... 100*sum(bddp_sim(:,phitry)./(1-rho*bdpdp_sim(:,phitry)) > bddp_data/(1-rho*phidp_data))/ntrial,... 100*sum(brdp_sim(:,phitry)./(1-rho*bdpdp_sim(:,phitry)) > (bddp_data+1-rho*phidp_data)/(1-rho*phidp_data))/ntrial,... 100*sum(bddp_sim(:,phitry)./(1-rho*bdpdp_sim(:,phitry)) > (brdp_data-1+rho*phidp_data)/(1-rho*phidp_data))/ntrial,... 100*sum(brdp_sim(:,phitry)./(1-rho*bdpdp_sim(:,phitry)) > brdp_data/(1-rho*phidp_data_impl))/ntrial,... 100*sum(bddp_sim(:,phitry)./(1-rho*bdpdp_sim(:,phitry)) > bddp_data/(1-rho*phidp_data_impl))/ntrial,... mean(brdp_sim(:,phitry)./(1-rho*bdpdp_sim(:,phitry))),mean(bddp_sim(:,phitry)./(1-rho*bdpdp_sim(:,phitry)))); fprintf(' %10s %10s %10s \n', 'phi', 'tr>data','td>data'); fprintf('%10.2f %10.2f %10.2f \n',... phiv(phitry),... 100*sum(trlr_sim(:,phitry) > trlr_data)/ntrial,... 100*sum(tdlr_sim(:,phitry) > tdlr_data)/ntrial); end; end; disp('******************'); disp('summary table'); % summary table fprintf('%10s & %10s & %10s & %10s & %10s & %10s & %10s & %10s \\\\ \n','phi','br','bdd','phi_br','blrmin','blrmax','s(dp)','1/2 life'); for phitry = 1:size(phiv,1); lrstats = [100*sum(brdp_sim(:,phitry)./(1-rho*bdpdp_sim(:,phitry)) > brdp_data/(1-rho*phidp_data))/ntrial; 100*sum(brdp_sim(:,phitry)./(1-rho*bdpdp_sim(:,phitry)) > (bddp_data+1-rho*phidp_data)/(1-rho*phidp_data))/ntrial; 100*sum(brdp_sim(:,phitry)./(1-rho*bdpdp_sim(:,phitry)) > brdp_data/(1-rho*phidp_data_impl))/ntrial]; fprintf('%10.3f & %10.2f & %10.2f & %10.2f & %10.2f & %10.2f & %10.2f & %10.1f \\\\ \n',... phiv(phitry),... 100*sum(brdp_sim(:,phitry)>brdp_data)/ntrial,... 100*sum(bddp_sim(:,phitry)>bddp_data)/ntrial,... 100*sum((bdpdp_sim(:,phitry)>phidp_data)&(brdp_sim(:,phitry)>brdp_data))/ntrial,... min(lrstats),max(lrstats),... (var_err_dp_data/(1-phiv(phitry)^2))^0.5,log(1/2)/log(phiv(phitry))); end; % two-way plots -- coefficients top = 1000; % how many data points to plot -- too many makes a fog marker_size = 8; line_width = 2; figure; for phitry = 1:size(phiv,1); subplot(2,2,2*(phitry-1)+1) plot(brdp_sim(1:top,phitry),bddp_sim(1:top,phitry),'.'); xlabel('b_r'); ylabel('b_d'); hold on; plot(brdp_data,bddp_data,'or','MarkerSize',marker_size,'MarkerFaceColor','r'); plot([brdp_data brdp_data]',[-0.4 0.3],'-r','LineWidth',line_width); plot([-0.2 0.6]',[bddp_data bddp_data]','-r','LineWidth',line_width); axis([-0.2 0.5 -0.35 0.2]); text(-0.1,0.15,[num2str(sum( (brdp_sim(:,phitry)bddp_data))/ntrial*100,'%5.1f') ' %']); text(0.25,0.15,[num2str(sum( (brdp_sim(:,phitry)>brdp_data)&(bddp_sim(:,phitry)>bddp_data))/ntrial*100,'%5.1f') ' %']); text(-0.1,-0.3,[num2str(sum( (brdp_sim(:,phitry)brdp_data)&(bddp_sim(:,phitry)tddp_data))/ntrial*100,'%5.1f') ' %']); text(2.5,2.2,[num2str(sum( (trdp_sim(:,phitry)>trdp_data)&(tddp_sim(:,phitry)>tddp_data))/ntrial*100,'%5.1f') ' %']); text(-1,-7,[num2str(sum( (trdp_sim(:,phitry)trdp_data)&(tddp_sim(:,phitry)phidp_data))/ntrial*100,'%5.1f') ' %']); text(0.15,1.00,[num2str(sum( (brdp_sim(:,phitry)>brdp_data)&(bdpdp_sim(:,phitry)>phidp_data))/ntrial*100,'%5.1f') ' %']); text(-0.05,0.75,[num2str(sum( (brdp_sim(:,phitry)brdp_data)&(bdpdp_sim(:,phitry)bddp_data)&(bdpdp_sim(:,phitry)bddp_data)&(bdpdp_sim(:,phitry)>phidp_data))/ntrial*100,'%5.1f') ' %']); text(1,-0.2,[num2str(sum( (bddp_sim(:,phitry)phidp_data))/ntrial*100,'%5.1f') ' %']); title(['b_d and \phi, \phi = ' num2str(phiv(phitry),'%5.2f')]); plot(phiv(phitry),-1+rho*phiv(phitry),'vg','MarkerSize',marker_size,'MarkerFaceColor','g'); bdline = (1-rho*phidp_data+bddp_data) + rho*[0.6 1.05]' -1; % bd: br = data, use inferred since (bd,phi) plot plot([0.6 1.05]',bdline,'--r','LineWidth',line_width); text(0.75,-0.2,'b_r','Color','r','BackgroundColor','w'); end; % new version with only br and phi figure; for phitry = 1:size(phiv,1); subplot(1,2,phitry) plot(brdp_sim(1:top,phitry),bdpdp_sim(1:top,phitry),'.'); xlabel('b_r'); ylabel('\phi'); hold on; plot(brdp_data,phidp_data,'or','MarkerSize',marker_size,'MarkerFaceColor','r'); plot([brdp_data brdp_data]',[0.6 1.05],'-r','LineWidth',line_width); plot([-0.1 0.4]',[phidp_data phidp_data]','-r','LineWidth',line_width); axis([-0.1 0.4 0.6 1.05]); text(-0.05,1.02,[num2str(sum( (brdp_sim(:,phitry)phidp_data))/ntrial*100,'%5.1f') ' %']); text(0.15,1.00,[num2str(sum( (brdp_sim(:,phitry)>brdp_data)&(bdpdp_sim(:,phitry)>phidp_data))/ntrial*100,'%5.1f') ' %']); text(-0.05,0.75,[num2str(sum( (brdp_sim(:,phitry)brdp_data)&(bdpdp_sim(:,phitry) lowlimit)&(brlr_sim(:,phitry) < hilimit); hist(brdp_sim(sel,phitry)./(1-rho*bdpdp_sim(sel,phitry)),50); h = findobj(gca,'Type','patch'); set(h,'FaceColor','y','EdgeColor','k') hold on; if ntrial == 50000; top = 3500; else; top = 280 end; plot(brdp_data/(1-rho*phidp_data)*[1 1]',[0 top]','-r','LineWidth',2); text(brdp_data/(1-rho*phidp_data),top*1.01,'Data '); set(gca,'Ytick',[]') % plot(1+bddp_data/(1-rho*phidp_data)*[1 1]',[0 220]','-r','LineWidth',3); % text(1+bddp_data/(1-rho*phidp_data),220,'b_d'); % f1 = 100*sum(brdp_sim(:,phitry)./(1-rho*bdpdp_sim(:,phitry)) > brdp_data/(1-rho*phidp_data))/ntrial; % f2 = 100*sum(bddp_sim(:,phitry)./(1-rho*bdpdp_sim(:,phitry)) > bddp_data/(1-rho*phidp_data))/ntrial; % text(1.3,280,[num2str(f1,'%5.2f') ' %']); % text(1.4,220,[num2str(f2,'%5.2f') ' %']); title(['\phi=' num2str(phi,'%5.2f')]); xlabel('b_r/(1-\rho\phi)'); %axis([lowlimit 2 0 400]); end; set(gcf,'PaperPosition',[0.25 2.5 8.0 3.0]); % default is [0.25 2.5 8.0 6.0] if print_graphs; print -depsc2 hist_longrun.eps; end; % t stats of the same figure; lowlimit = -inf; % cuts off a very few big negatives that spoil histogram for phitry = 1:size(phiv,1); phi = phiv(phitry); subplot(1,2,phitry) sel = trlr_sim(:,phitry) > lowlimit; hist(trlr_sim(sel,phitry),50); h = findobj(gca,'Type','patch'); set(h,'FaceColor','y','EdgeColor','k') hold on; plot(trlr_data*[1 1]',[0 280]','-r','LineWidth',2); % text(trlr_data,280,'b_r '); %f1 = 100*sum(trlr_sim(:,phitry) > trlr_data)/ntrial; %text(1.3,280,[num2str(f1,'%5.2f') ' %']); title(['\phi=' num2str(phi,'%5.2f')]); xlabel('t r long run'); %axis([lowlimit 2 0 400]); end; set(gcf,'PaperPosition',[0.25 2.5 8.0 3.0]); % default is [0.25 2.5 8.0 6.0] if print_graphs; print -depsc2 hist_longrun_t.eps; end; % -------------------------- % comparison graph in br-phi space including likelihood % ----------------------- phi_plot = (0.6:0.002:1)'; %k = (0:0.025:0.2); k = (exp(chi2inv([0.95],1)/76)-1) * std(err_r)^2/std(err_dp)^2 ; % add [0.90 0.95 0.99] to see more p values phi_plot = phi_plot*ones(1,size(k,2)); k = ones(size(phi_plot,1),1)*k; br_plot = (k.*(1-phi_plot.^2)).^0.5; br_plot_neg = -(k.*(1-phi_plot.^2)).^0.5; figure; % asymptotic LR regions plot(br_plot,phi_plot,'-m','LineWidth',2); hold on plot(br_plot_neg,phi_plot,'-m','LineWidth',2); % legend of % prob values for LR lines phi_plot = 0.75; k = (exp(chi2inv([0.95 ],1)/76)-1) * std(err_r)^2/std(err_dp)^2 ; br_plot = (k.*(1-phi_plot.^2)).^0.5; text(br_plot',phi_plot,['LR 5%'],'BackgroundColor','w','HorizontalAlignment','center','color','m') % simulated data plot(brdp_sim(1:top,1),bdpdp_sim(1:top,1),'.','MarkerSize',3); % br = data line plot([brdp_data brdp_data]',[0.6 1.05],'-r','LineWidth',line_width); text(0.09,0.7,'b_r','Color','r','BackgroundColor','w'); % bd = data line philine1 = 1/rho*(- [-0.1 0.4]' + 1 + (brdp_data - 1 + rho*phidp_data)); % phi = (-br+ 1 + bd)/rho since (br, phi) graph use bd implied by % br phi rather than exact, then it will go through points. plot([-0.1 0.4]',philine1,'--r','LineWidth',line_width); % bd event, keeping bddp at bddpdata text(0.3,0.70,'b_d','Color','r','BackgroundColor','w'); % data point and null plot(brdp_data,phidp_data,'or','MarkerSize',marker_size,'MarkerFaceColor','r'); plot(0,phiv(1),'vg','MarkerSize',marker_size,'MarkerFaceColor','g'); % long run regressions phi_plot1 = [0.5:0.01:1.1]'; br_plot_last = brdp_data.*(1-rho*phi_plot1)./(1-rho*phidp_data); plot(br_plot_last,phi_plot1,'-b','LineWidth',line_width) text(0.32,0.75,'b_r^{lr}','Color','b','BackgroundColor','w'); axis([-0.1 0.37 0.6 1.03]); xlabel('b_r') ylabel('\phi'); set(gca,'Xtick',[-0.1:0.1:0.4]); set(gca,'Ytick',[0.6:0.1:1]); print -depsc2 summary_region.eps % -------- % distributions of direct long run estimates. % ---------------------- figure; subplot(2,2,1); hist(lr5(:,1),50); hold on; plot(lr5_data*[1 1]',[0 330]','-r','LineWidth',3); text(lr5_data+0.3,200,[num2str(sum(lr5(:,1)>lr5_data)/size(lr5(:,1),1)*100,'%3.1f') '%']); title('5 year returns '); axis([-0.5 1.5 0 350]); set(gca,'Ytick',[]) xlabel('Regression Coefficient') subplot(2,2,2); hist(lr10(:,1),50); hold on; plot(lr10_data*[1 1]',[0 330]','-r','LineWidth',3); text(lr10_data+0.3,200,[num2str(sum(lr10(:,1)>lr10_data)/size(lr10(:,1),1)*100,'%3.1f') '%']); title('10 year returns'); axis([-1.5 2 0 350]); set(gca,'Ytick',[]) xlabel('Regression Coefficient') subplot(2,2,3); hist(lr15(:,1),50); hold on; plot(lr15_data*[1 1]',[0 330]','-r','LineWidth',3); text(lr15_data+0.3,200,[num2str(sum(lr15(:,1)>lr15_data)/size(lr15(:,1),1)*100,'%3.1f') '%']); title('15 year returns '); axis([-1.5 2.5 0 350]); set(gca,'Ytick',[]) xlabel('Regression Coefficient') subplot(2,2,4); hist(lr20(:,1),50); hold on; plot(lr20_data*[1 1]',[0 330]','-r','LineWidth',3); text(lr20_data+0.3,200,[num2str(sum(lr20(:,1)>lr20_data)/size(lr20(:,1),1)*100,'%3.1f') '%']); title('20 year returns'); axis([-1.3 2.5 0 350]); set(gca,'Ytick',[]) xlabel('Regression Coefficient') print -depsc2 longrun_direct.eps disp('direct long run regressions, finite horizons'); disp('percentage of observations with coefficient greater than data'); disp(' weighted long-run sum(rho^j x_t+j)') fprintf(' phi = %10s %10.2f %10.2f \n',' ', phiv(1), phiv(2)); fprintf(' 1 year returns %10.2f %10.2f %10.2f \n',lr1_data,sum(lr1>lr1_data)/size(lr1,1)*100); fprintf(' 5 year returns %10.2f %10.2f %10.2f \n',lr5_data,sum(lr5>lr5_data)/size(lr5,1)*100); fprintf('10 year returns %10.2f %10.2f %10.2f \n',lr10_data,sum(lr10>lr10_data)/size(lr10,1)*100); fprintf('15 year returns %10.2f %10.2f %10.2f \n',lr15_data,sum(lr15>lr15_data)/size(lr15,1)*100); fprintf('20 year returns %10.2f %10.2f %10.2f \n',lr20_data,sum(lr20>lr20_data)/size(lr20,1)*100);% fprintf(' 1 year dividends %10.2f %10.2f %10.2f \n',ld1_data,sum(ld1>ld1_data)/size(ld1,1)*100); fprintf(' 5 year dividends %10.2f %10.2f %10.2f \n',ld5_data,sum(ld5>ld5_data)/size(ld5,1)*100); fprintf('10 year dividends %10.2f %10.2f %10.2f \n',ld10_data,sum(ld10>ld10_data)/size(ld10,1)*100); fprintf('15 year dividends %10.2f %10.2f %10.2f \n',ld15_data,sum(ld15>ld15_data)/size(ld15,1)*100); fprintf('20 year dividends %10.2f %10.2f %10.2f \n',ld20_data,sum(ld20>ld20_data)/size(ld20,1)*100); disp(' unweighted long-run sum( x_t+j)') fprintf(' phi = %10.s %10.2f %10.2f \n',' ',[phiv(1) phiv(2)]); fprintf(' 1 year returns %10.2f %10.2f %10.2f \n',lru1_data,sum(lru1>lru1_data)/size(lru1,1)*100); fprintf(' 5 year returns %10.2f %10.2f %10.2f \n',lru5_data,sum(lru5>lru5_data)/size(lru5,1)*100); fprintf('10 year returns %10.2f %10.2f %10.2f \n',lru10_data,sum(lru10>lru10_data)/size(lru10,1)*100); fprintf('15 year returns %10.2f %10.2f %10.2f \n',lru15_data,sum(lru15>lru15_data)/size(lru15,1)*100); fprintf('20 year returns %10.2f %10.2f %10.2f \n',lru20_data,sum(lru20>lru20_data)/size(lru20,1)*100); fprintf(' 1 year dividends %10.2f %10.2f %10.2f \n',ldu1_data,sum(ldu1>ldu1_data)/size(ldu1,1)*100); fprintf(' 5 year dividends %10.2f %10.2f %10.2f \n',ldu5_data,sum(ldu5>ldu5_data)/size(ldu5,1)*100); fprintf('10 year dividends %10.2f %10.2f %10.2f \n',ldu10_data,sum(ldu10>ldu10_data)/size(ldu10,1)*100); fprintf('15 year dividends %10.2f %10.2f %10.2f \n',ldu15_data,sum(ldu15>ldu15_data)/size(ldu15,1)*100); fprintf('20 year dividends %10.2f %10.2f %10.2f \n',ldu20_data,sum(ldu20>ldu20_data)/size(ldu20,1)*100); disp('long run regressions implied by short run coefficients'); disp(' weighted long-run sum(rho^j x_t+j)') fprintf(' phi = %10s %10.2f %10.2f \n',' ', phiv(1), phiv(2)); fprintf(' 1 year returns %10.2f %10.2f %10.2f \n',brdp_data,... sum(brdp_sim>brdp_data)/size(brdp_sim,1)*100); fprintf(' 2 year returns %10.2f %10.2f %10.2f \n',brdp_data*(1+rho*phidp_data),... sum(brdp_sim.*(1+rho*bdpdp_sim)>brdp_data*(1+rho*phidp_data))/size(brdp_sim,1)*100); fprintf(' 5 year returns %10.2f %10.2f %10.2f \n',brdp_data*(1-rho^5*phidp_data^5)/(1-rho*phidp_data),... sum(brdp_sim.*(1-rho^5*bdpdp_sim.^5)./(1-rho*bdpdp_sim)>... brdp_data*(1-rho^5*phidp_data^5)/(1-rho*phidp_data))/size(brdp_sim,1)*100); fprintf('10 year returns %10.2f %10.2f %10.2f \n',brdp_data*(1-rho^10*phidp_data^10)/(1-rho*phidp_data),... sum(brdp_sim.*(1-rho^10*bdpdp_sim.^10)./(1-rho*bdpdp_sim)>... brdp_data*(1-rho^10*phidp_data^10)/(1-rho*phidp_data))/size(brdp_sim,1)*100); fprintf('15 year returns %10.2f %10.2f %10.2f \n',brdp_data*(1-rho^15*phidp_data^15)/(1-rho*phidp_data),... sum(brdp_sim.*(1-rho^15*bdpdp_sim.^15)./(1-rho*bdpdp_sim)>... brdp_data*(1-rho^15*phidp_data^15)/(1-rho*phidp_data))/size(brdp_sim,1)*100); fprintf('20 year returns %10.2f %10.2f %10.2f \n',brdp_data*(1-rho^20*phidp_data^20)/(1-rho*phidp_data),... sum(brdp_sim.*(1-rho^20*bdpdp_sim.^20)./(1-rho*bdpdp_sim)>... brdp_data*(1-rho^20*phidp_data^20)/(1-rho*phidp_data))/size(brdp_sim,1)*100); fprintf('inf year returns %10.2f %10.2f %10.2f \n',brdp_data/(1-rho*phidp_data),... sum(brdp_sim./(1-rho*bdpdp_sim)>... brdp_data/(1-rho*phidp_data))/size(brdp_sim,1)*100); disp(' unweighted long-run sum(x_t+j)') fprintf(' phi = %10s %10.2f %10.2f \n',' ', phiv(1), phiv(2)); fprintf(' 1 year returns %10.2f %10.2f %10.2f \n',brdp_data,... sum(brdp_sim>brdp_data)/size(brdp_sim,1)*100); fprintf(' 2 year returns %10.2f %10.2f %10.2f \n',brdp_data*(1+phidp_data),... sum(brdp_sim.*(1+bdpdp_sim)>brdp_data*(1+phidp_data))/size(brdp_sim,1)*100); fprintf(' 5 year returns %10.2f %10.2f %10.2f \n',brdp_data*(1-phidp_data^5)/(1-phidp_data),... sum(brdp_sim.*(1-bdpdp_sim.^5)./(1-bdpdp_sim)>... brdp_data*(1-phidp_data^5)/(1-phidp_data))/size(brdp_sim,1)*100); fprintf('10 year returns %10.2f %10.2f %10.2f \n',brdp_data*(1-rho^10*phidp_data^10)/(1-phidp_data),... sum(brdp_sim.*(1-bdpdp_sim.^10)./(1-bdpdp_sim)>... brdp_data*(1-phidp_data^10)/(1-phidp_data))/size(brdp_sim,1)*100); fprintf('15 year returns %10.2f %10.2f %10.2f \n',brdp_data*(1-rho^15*phidp_data^15)/(1-phidp_data),... sum(brdp_sim.*(1-bdpdp_sim.^15)./(1-bdpdp_sim)>... brdp_data*(1-phidp_data^15)/(1-phidp_data))/size(brdp_sim,1)*100); fprintf('20 year returns %10.2f %10.2f %10.2f \n',brdp_data*(1-rho^20*phidp_data^20)/(1-phidp_data),... sum(brdp_sim.*(1-bdpdp_sim.^20)./(1-bdpdp_sim)>... brdp_data*(1-phidp_data^20)/(1-phidp_data))/size(brdp_sim,1)*100); fprintf('inf year returns %10.2f %10.2f %10.2f \n',brdp_data/(1-phidp_data),... sum(brdp_sim./(1-bdpdp_sim)>... brdp_data/(1-phidp_data))/size(brdp_sim,1)*100); % br, phi plot with lines to explain long run estimates figure; phitry = 1; plot(brdp_sim(1:top,phitry),bdpdp_sim(1:top,phitry),'.'); xlabel('b_r'); ylabel('\phi'); hold on; plot(brdp_data,phidp_data,'or','MarkerSize',marker_size,'MarkerFaceColor','r'); axis([-0.05 0.4 0.65 1.05]); title('b_r and \phi, with long run regressions ' ); plot(0,phiv(phitry),'vg','MarkerSize',marker_size,'MarkerFaceColor','g'); k1 = [ 1 5 10 20]'; Nk = size(k1,1); phi_plot1 = [0.5:0.01:1.1]'; Nphi = size(phi_plot1,1); k = ones(Nphi,1)*k1'; phi_plot = phi_plot1*ones(1,Nk); br_plot = brdp_data.*(1-rho.^k.*phidp_data.^k)./(1-rho.^k.*phi_plot.^k).*(1-rho*phi_plot)./(1-rho*phidp_data); plot(br_plot,phi_plot,'-r','LineWidth',line_width); % infinity as a special case br_plot_last = brdp_data.*(1-rho*phi_plot1)./(1-rho*phidp_data); plot(br_plot_last,phi_plot1,'-r','LineWidth',line_width) br_plot_unw = brdp_data.*(1-phi_plot1)./(1-phidp_data); plot(br_plot_unw,phi_plot1,'--r','LineWidth',line_width) phi_plot_2 = 0.7; br_plot2 = brdp_data.*(1-rho.^k1.*phidp_data.^k1)./(1-rho.^k1.*phi_plot_2.^k1).*(1-rho*phi_plot_2)./(1-rho*phidp_data); text(br_plot2,phi_plot_2*ones(Nk,1),[' 1';' 5';'10';'20'],'BackgroundColor',[1 1 1],'HorizontalAlignment','center'); br_plot_last2 = brdp_data.*(1-rho*phi_plot_2)./(1-rho*phidp_data); text(br_plot_last2,phi_plot_2,'\infty','BackgroundColor',[1 1 1],'HorizontalAlignment','center'); phi_plot_3 = 0.82; br_plot_unw3 = brdp_data.*(1-phi_plot_3)./(1-phidp_data); text(br_plot_unw3+0.03,phi_plot_3,'\infty, unweighted','BackgroundColor',[1 1 1],'HorizontalAlignment','center'); set(gca,'xtick',[0:0.1:0.4]); set(gca,'ytick',[0.7:0.1:1]); figheight = 4; figwidth = 6; set(gcf, 'PaperPositionMode', 'manual'); set(gcf, 'PaperUnits', 'inches'); set(gcf, 'PaperPosition', [0.25 2.5 figwidth figheight]); if print_graphs print -depsc2 explain_long_run.eps end; % same for unweighted. Isn't different enough to show, though the infinity case is cool since the line % is "too steep". figure; plot(brdp_sim(1:top,phitry),bdpdp_sim(1:top,phitry),'.'); xlabel('b_r'); ylabel('\phi'); hold on; plot(brdp_data,phidp_data,'or','MarkerSize',marker_size,'MarkerFaceColor','r'); axis([-0.05 0.4 0.65 1.05]); title('b_r and \phi, with unweighted long run regressions ' ); plot(0,phiv(phitry),'vg','MarkerSize',marker_size,'MarkerFaceColor','g'); k1 = [ 1 2 5 10 15 20]'; Nk = size(k1,1); phi_plot1 = [0.5:0.01:1.1]'; Nphi = size(phi_plot1,1); k = ones(Nphi,1)*k1'; phi_plot = phi_plot1*ones(1,Nk); br_plot = brdp_data.*(1-phidp_data.^k)./(1-phi_plot.^k).*(1-phi_plot)./(1-phidp_data); plot(br_plot,phi_plot,'-r','LineWidth',line_width); % infinity as a special case br_plot_last = brdp_data.*(1-phi_plot1)./(1-phidp_data); plot(br_plot_last,phi_plot1,'-r','LineWidth',line_width) phi_plot_2 = 0.7; br_plot2 = brdp_data.*(1-phidp_data.^k1)./(1-phi_plot_2.^k1).*(1-phi_plot_2)./(1-phidp_data); text(br_plot2,phi_plot_2*ones(Nk,1),[' 1';' 2';' 5';'10';'15';'20'],'BackgroundColor',[1 1 1],'HorizontalAlignment','center'); br_plot_last2 = brdp_data.*(1-phi_plot_2)./(1-phidp_data); text(br_plot_last2,phi_plot_2,'\infty','BackgroundColor',[1 1 1],'HorizontalAlignment','center'); % not for paper, but interesting: joint dist of t stat and of estimate along with phi % makes the point that the estimate is controlling for phi figure subplot(1,2,1); plot(trlr_sim(1:top,1),bdpdp_sim(1:top,1),'.'); xlabel('t long run'); ylabel('\phi'); subplot(1,2,2) plot(brlr_sim(1:top,1),bdpdp_sim(1:top,1),'.'); xlabel('br long run'); ylabel('phi') % find ``confidence ellipses'' if 0; % takes a while, and not used in paper. for phitry = 1:size(phiv,1); step = 0.02; brgrid = (-0.2:step:0.5)'; brgrid = sort([brdp_data;brgrid]); bdgrid = (-0.4:step:0.15)'; bdgrid = sort([bddp_data;bdgrid]); jointdist = 0*bdgrid*brgrid'; histdist = 0*bdgrid*brgrid'; s = 0.01; for i = 1:size(bdgrid,1); for j = 1:size(brgrid,1); jointdist(i,j) = sum(exp( -(bddp_sim(:,phitry)-bdgrid(i)).^2/(2*s^2) - (brdp_sim(:,phitry)-brgrid(j)).^2/(2*s^2))); %kernel histdist(i,j) = sum( (abs(bddp_sim(:,phitry)-bdgrid(i)) < step/2-1e-6) & (abs(brdp_sim(:,phitry)-brgrid(j)) < step/2-1e-6)); %kernel end; end; jointdist = jointdist/sum(sum(jointdist)); % find level of jointdist such that 5% of joint distrbution is in boxes below this level % use regular histogram for the 5%; jointdist being a kernel spreads out too wide figure lev = findlev(histdist,jointdist,5); disp('size'); disp(100*sum(sum(histdist(jointdist lev; % a mesa for the included region contour(brgrid,bdgrid,region,1,'LineWidth',3); hold on; plot(brdp_sim(:,phitry),bddp_sim(:,phitry),'.') plot(brdp_data,bddp_data,'or','MarkerSize',5,'MarkerFaceColor','r'); rindx = find(brgrid==brdp_data); dindx = find(bdgrid==bddp_data); condist_r = jointdist(dindx,:)'/sum(jointdist(dindx,:)'); condist_d = jointdist(:,rindx)/sum(jointdist(:,rindx)); end; % plot conditional distribution of r, d given the other. % only94 = 1; % do only phi = 0.94 results maxphi = size(phiv,1); if only94; maxphi = 1; end; figure; for phitry = 1:maxphi; brgrid = (-0.1:0.01:0.4)'; brgrid = sort([brdp_data;brgrid]); bdgrid = (-0.3:0.01:0.1)'; bdgrid = sort([bddp_data;bdgrid]); jointdist = bdgrid*brgrid'; s = 0.02; for i = 1:size(bdgrid,1); for j = 1:size(brgrid,1); jointdist(i,j) = sum(exp( -(bddp_sim(:,phitry)-bdgrid(i)).^2/(2*s^2) - (brdp_sim(:,phitry)-brgrid(j)).^2/(2*s^2))); end; end; jointdist = jointdist/sum(sum(jointdist)); rindx = find(brgrid==brdp_data); dindx = find(bdgrid==bddp_data); condist_r = jointdist(dindx,:)'/sum(jointdist(dindx,:)'); condist_d = jointdist(:,rindx)/sum(jointdist(:,rindx)); if only94 == 0; if phitry == 1; subplot(2,2,1); elseif phitry == 2; subplot(2,2,3); end; elseif only94 == 1; subplot(1,2,1); end; plot(brgrid,condist_r,'LineWidth',2); hold on; plot(brdp_data*[1 1]', [0 .08], '-r','LineWidth',3); text(brdp_data+0.005,0.08,'b_r in data'); text(0,0.05,[int2str(sum(condist_r(1:rindx))/sum(condist_r)*100) ' %']); text(0.25,0.05,[int2str(sum(condist_r(rindx+1:end))/sum(condist_r)*100) ' %']); xlabel('b_r'); axis([-0.1 0.4 0 0.1]); set(gca,'Ytick',[]); if only94 == 1; title('f( b_r | b_d )'); else; title(['f( b_r | b_d ), \phi = ' num2str(phiv(phitry),'%5.2f')]); end; if only94 == 0; if phitry == 1; subplot(2,2,2); elseif phitry == 2; subplot(2,2,4); end; elseif only94 == 1; subplot(1,2,2); end; plot(bdgrid,condist_d,'LineWidth',2); hold on; plot(bddp_data*[1 1]',[0 0.10],'-r','LineWidth',3); text(bddp_data+0.005,0.10,'b_d in data'); xlabel('b_d') text(-0.2,0.05,[int2str(sum(condist_d(1:dindx))/sum(condist_d)*100) ' %']); text(0.04,0.05,[int2str(sum(condist_d(dindx+1:end))/sum(condist_d)*100) ' %']); axis([-0.3 0.1 0 0.12]); set(gca,'Ytick',[]); if only94 == 0; title(['f( b_d | b_r ), \phi = ' num2str(phiv(phitry),'%5.2f')]); end; if only94 == 1; title('f( b_d | b_r )'); end; end; if only94 == 1; set(gcf,'PaperPosition',[0.25 2.5 8.0 3.0]); % default is [0.25 2.5 8.0 6.0] end; if print_graphs; print -depsc2 conditional_dist_b.eps end; end; % Goyal welch if doGW; lowlimit = -2.5; uplimit = 1.5; % limts for graphing histogram figure; subplot(1,2,1) sel = (msediff(:,1)>lowlimit)&(msediff(:,1)lowlimit)&(msediff(:,1)