function [m,s]=estimateMultModelParams(stockTimeSeries,indexTimeSeries,LookbackPeriod,indexReturn,riskFreeRate)

swr=stockTimeSeries(1:LookbackPeriod)./stockTimeSeries(2:LookbackPeriod+1)-1;

iwr=indexTimeSeries(1:LookbackPeriod)./indexTimeSeries(2:LookbackPeriod+1)-1;

covmatx=cov(swr,iwr);

beta=covmatx(1,2)/covmatx(2,2);

fcastR=riskFreeRate/LookbackPeriod+beta*(indexReturn/LookbackPeriod-riskFreeRate/LookbackPeriod);

s=std(log(stockTimeSeries(1:LookbackPeriod)./stockTimeSeries(2:LookbackPeriod+1)));

m=log(1+fcastR)-s^2/2;

end

If you define and input data into those variables, you can use it to estimate the probability density of a function of  a variable subjected to normal distribution, with a mean of m, and a standard deviation of s. Here is an example:

>>x=-4:.1:4;

>>m=.06;

>>s=.02;

>>T=5;

>>By=1./(1+(m+s*x)/2).^(2*T);

>>Plot(By, (1/(s*T))*normpdf(x)./(By.^(1+1/(2*T)))), title(‘bond price distribution with a yield subject to normal distribution.’)