clear all
V=1;
N=10000;

while V<=N

T=390;   % Number of time intervals (390 is the number of minutes in a trading day)
  
Sign=round(rand(1,T));   % Vector that determines the direction of price change
S=-1*(Sign<1);
Sign=Sign+S;

Delta=round(rand(1,T));   % Random price change vector (e.g. [0 -1 -1 0 1 0 1])
Delta=Delta.*Sign;

m=1;
P(1)=40;
P(T)=40;
t(1)=1;
VolatilityFunction(1)=20;
PriceChangeIncrement=0.01;   % When the price changes it will change by this amount (Equals damping term/mass)
Freq(1)=0.4;   % Initial frequency of volatility oscillation
PriceChange(1)=0;
MovingAvgPriceChange(1)=0;

for m=1:T
Freq(m+1)=Freq(1)*exp(-m/T);   % Decreases oscillation frequency as time moves forward
VF=VolatilityFunction(1)+VolatilityFunction(m)*(P(T)/P(1))*exp(-PriceChangeIncrement*m)*sin(Freq(m+1)*m);   % PT/P1 can be thought of as "springiness"
   VolatilityFunction(m+1)= max(VF,-VF);
   P(m+1)=P(m)+PriceChangeIncrement*Delta(m)*VolatilityFunction(m)+(P(T)-P(1))/T;
   PriceChange(m+1)=(P(m+1)-P(m))/P(m);
   MovingAvgPriceChange(m+1)=mean(PriceChange(max(m-round(2*pi/Freq(m)),1):m));   % Moving average time period changes with frequency of oscillation 
   t(m+1)=m;
   m=m+1;
end

for q=1:T
PricePoints(V,q)=P(q);
end

V=V+1;
end


for r=1:T
MeanPrice(r)=mean(PricePoints(1:N,r));
MedianPrice(r)=median(PricePoints(1:N,r));
PriceStdDev(r)=std(PricePoints(1:N,r));
end

MeanVolatility=mean(VolatilityFunction)
MeanVolatility=mean(VolatilityFunction)*ones(1,T+1);

subplot(4,1,1), hist(PricePoints(:,T),50), title('Sample Distribution of Closing Prices')

subplot(4,1,2), plot(t,P,'k'), AXIS([0 T (min(PricePoints(N,1:T)-PriceStdDev(1:T))-5) (max(PricePoints(N,1:T)+PriceStdDev(1:T))+5)]), title('Single Realization of Asset Random Walk')
hold on
subplot(4,1,2), plot(0:T-1,(PricePoints(N,1:T)+PriceStdDev(1:T)),'r:')
subplot(4,1,2), plot(0:T-1,(PricePoints(N,1:T)-PriceStdDev(1:T)),'r:')

subplot(4,1,3), plot(t,PriceChange,'g'), AXIS([0 T (min(PriceChange)-0.005) (max(PriceChange)+0.005)]), title('Single Realization of Percentage Change in Asset Price')
hold on
subplot(4,1,3), plot(t,MovingAvgPriceChange,'k')

subplot(4,1,4), plot(t,VolatilityFunction,'r'), , AXIS([0 T 0 (max(VolatilityFunction)+5)]), title('Single Realization of Asset Price Volatility')
hold on
subplot(4,1,4), plot(t,MeanVolatility,'k:')

figure
plot(t,Freq)
grid on
title('Frequency of Oscillation in Volatility vs. Time')
xlabel('Time')
ylabel('Frequency of Oscillation')