Partition function in Matlab - is there something missing in my code?

Can anyone see if there is something wrong in my matlab code? My objective is to replicate this formula: http://i.stack.imgur.com/JANLf.jpg q can take value 1,2,3 and 5. I constructed my vector Xt where each element are a cumulative sum of log(1+return) at each time (t) - for stock returns - first element is normalized to log(1).

Then to compute each element Sq(T,delta t) for the four values of q this is my matlab code:

for j=1:length(dt);
E=Xt(1:dt(j):end);
EE=diff(E(2:end));
EEE=diff(E(1:end-1));
Sqone(j)=sum(abs(EE-EEE).^1);
Sqtwo(j)=sum(abs(EE-EEE).^2);
Sqthree(j)=sum(abs(EE-EEE).^3);
Sqfive(j)=sum(abs(EE-EEE).^5); end;
Is there something wrong in the code above? I am asking this because I know there is something wrong since I am not getting the expected results. I am convinced that it is due to my code posted above.

the vector dt is a vector that goes from 1 to high number - depending on the size of Xt. But my vector dt is not the problem.

Thank you for all your help!

Answer 1

[Andrei Bobrov]

Sq = @(x,q)sum(abs(diff(x)).^q)

eg

Sq(Xt,5)

EDIT

Sq = @(x,q)sum(bsxfun(@power,abs(diff(x(:))),q(:)'));

eg

Sq(Xt,[1:3 5])

EDIT2

eg

Xt = randi(15,15,1);
dt = 1:14;
q = [1:3,5];
out = zeros(numel(dt),numel(q));
for j1 = 1:numel(dt)
     out(j1,:) = sum(bsxfun(@power,abs(diff(Xt(1:dt(j1):end))),q),1);
end

variant without loop for...end

out2 = cell2mat(cellfun(@(x) sum(bsxfun(@power,abs(diff(x(:))),q),1),arrayfun(@(x)Xt(1:x:end),dt(:),'un',0),'un',0));

[Charles Martineau]

Andrei thanks for this superb formula... it improves my code but then I still can't solve my problem. I posted my entire code below for a Brownian motion while including your code. I don't know if you can spot the error! THANKS

Answer 2

[Charles Martineau]

Hi Andrei, I don't really understand. Is this code part of a loop? do I have to implement values like for instance; Sq = @(x,q)sum(abs(diff(x)).^q) becomes Sq = @(Xt,5)sum(abs(diff(Xt)).^5)? Thanks for the help

[Andrei Bobrov]

No, this is a complete solution. 'Sq' - handle function

[Charles Martineau]

Ok now I see that this formula gives me one number or one point. But I need to construct a Sq vector 1 X size(dt).

[Andrei Bobrov]

What is it 'dt' and 'Xt'. Please give simple numeric example.

[Charles Martineau]

dt (delta t) is a vector that takes different values for the change in t. So for instance dt = [1 2 3 4 5 6 7 8.... ] size 1 X length(Xt). Xt is a vector that takes the log price series lnP(t) - lnP(0) - size 1 X size(Price vector).

Thanks!

Answer 3

[Charles Martineau]

Andrei thanks for this superb formula! It improves my code but for some reason I still generate the same or similar output has my previous formula.

If I use your code above in a simulation of a Brownian motion for the vector Xt of log(price)

T=1000; 
Drift_annual=.10; %pick a number
Volatility_annual=.40; %pick a number
Drift_dayly=Drift_annual/T;
Volatility_dayly=Volatility_annual/sqrt(T);
Drift_mean=Drift_dayly-(.5*(Volatility_dayly^2));
Normal=randn(1,T);
Pvec=1; %Price vector

for t=1:T;
    LogR(t)=Drift_mean+Normal(t)*Volatility_dayly;
end;

for t=2:T;
Pvec(t)=Pvec(t-1)*exp(LogR(t-1));
end;

lnp=log(Pvec);
Xt=1;
for m=2:length(Pvec);
    Xt(m)=lnp(m)-lnp(1);
end;
Xt=Xt'

N=max(size(Xt));
minobs=30; %minimum observations

%Define interval dt%
dlogdt=.05;
dt=round(exp([0:dlogdt:log(N/minobs)]));
dt=dt([1, 1+find(dt(1:length(dt)-1)~=dt(2:length(dt)))]);
ddt=1:length(dt);
qvec=[1:3 5];
nq=length(qvec);
ndt=length(dt);
cut1=1;
cut2=ndt;

out = zeros(numel(dt),numel(qvec));
for j1 = 1:numel(dt)
     out(j1,:) = sum(bsxfun(@power,abs(diff(Xt(1:dt(j1):end))),qvec),1);
end

partitionfn=(out')./(ones(nq,1)*dt); %Here I make sure that I rescale the values to start at 10^0 when plotted on a loglog graph%
pfngrap=(partitionfn./(partitionfn(:,cut1)*ones(1,ndt)));

figure(1);
loglog(dt(cut1:cut2),pfngrap(:,cut1:cut2));

hold on;
xmin=dt(cut1);
xmax=max(dt);
brownmat=exp((((qvec/2)-1)'*([0 log(xmax/xmin)])));
plot([xmin xmax],brownmat,'r:');
hold off;
ymin=min(brownmat(:,2));
ymax=max(brownmat(:,2));
axis([xmin xmax ymin ymax]);

my generated "out" vector elements, when plotted on a loglog should fall on the red lines in the generated graph but they don't.... . This is the real way for me to make sure that I wrote the code correctly. If you simply copy paste the code above and run it you will see that there is something wrong. Can you spot it? I have been losing many days trying to figure it out... Thanks again for all the much appreciated help!