s-distributionprograms

Computer programs For S-Distributions
      Workshop Main Page
a. Calculating the fixing parameter "τ(m,ρ)" ( in the sense that
it fixes the mixture to be density).
Search The body function of the root and suitable r , as it
explained in the paper, the program are in the file named as
“Tao.txt” and “FoundR”. The main program is in the file”
Positive Root” which contain these parts:
Program Environment: MATLAB 7.8.0(R2009a)
Inputs:
“N, P”, which are “m, ρ”, the mixture order parameter and the
right side skewness parameter, respectively in the paper.
Outputs:
Rmin: which is the fixing parameter "τ(m,ρ)" .

a1. For the positive Root:

function [ Rmin ]= Positive_Root(N,p,rstep)

L=[];
for i= 1:N+1
   l=(0:1/i:1);
   l=l(1:end);
   L=[L,[l;(0:i);zeros(size(l))+i-1]];
end

index_p=1;
index_k=2;
index_N=3;

[b m n]=unique(L(index_p,:));
[a in]=sort(L(index_p,:));

QL=L(:,m);
NQ=L(:,in);
PNQ=NQ(index_p,:)>0;
NQL=NQ(:,PNQ);
L_in=find(QL(index_p,:)<p);

is=0;
[X R tao]=FoundR(NQL,QL,p,N,rstep,is);

Xmin=X;
Rmin=R;
MXmin=[X];
MRmin=[R];
for is=1:size(L_in,2)-2
  [x r]=FoundR(NQL,QL,p,N,rstep,is);
  MXmin=[x MXmin];
  MRmin=[r MRmin];
  if r<Rmin
      Rmin=r;
  end
end
%
end

a2. Tao:
function tao=Tao(n,k,p)
Rn=inline('sin((k.*pi)/(n+1)) ./ sin(p.*pi - (k.*pi)/(n+1))');
tao=Rn(k,n,p);
end

a3. FoundR:

function [X R tao]=FoundR(NQL,QL,p,N,rstep,is)

   index_p=1;
   index_k=2;
   index_N=3;

   L_in=find(QL(index_p,:)<p);
   lp=QL(index_p,L_in(end));

   U_in=find(QL(index_p,:)>p);
   if is==0
       U_in=find(NQL(index_p,:)==lp);%QL(index_p,U_in(1)));
   else
       U_in=find(NQL(index_p,:)==QL(index_p,L_in(end-is)));
   end
   up=NQL(index_p,U_in(end)+1);
   Up=NQL(:,U_in(end));

   Q=NQL(:,find(NQL(index_p,:)<up));

   T=[];
   for n=1:N
       if mod(sum(Q(index_N,:)==n),2)==1
           T=[T n];
       end
   end

   n=Up(index_N);
   k=Up(index_k);

   r=rstep:rstep:10;
   tao=repmat(Tao(n,k,p),[size(r)]);
   rtao=r./tao;

   SRT=zeros(size(r));
   for i=1:size(T,2)
       SRT=SRT+rtao.^T(i);
   end

   SRT=SRT.*(SRT<1);
   [srt indexr]=max(SRT);
   R=r(indexr);
   X=rtao(indexr);
   tao=tao(1);
end

b. MLE program for parameters estimation:

Programs Environment: MATLAB 7.8.0(R2009a)

For calculating MLE's parameters, firstly we have to define the
body function of the Log likelihood for the S distribution,
which have three parts in the files named as (f1.txt, avg_f1.txt,
f2.txt) respectively. Then in the program named as, main7.txt ,
we estimate the parameter by using MLE
method, this program has these parts:
Program Environment: MATLAB 7.8.0(R2009a)
Inputs:
finp1='tr.xls'; which is the matrix of "m " and " ρ" and the
fixing parameter "τ(m,ρ)" ( in the sense that it fixes the mixture
to be density).
finp2='AMEX.xls'; the file of data that we can used
dat: for reading data from excel file
m: which is the mixture order parameter
Outputs:
The MLE's estimated parameters, "rh, ah, bh, rhoh " which are
" r, α, λ, ρ" respectively, in the paper.

b1. f1.txt:

function y=f1(t,n,rho,a,b)
y=0;
for j=1:(n+1)
y=y+nchoosek(n+1,j).*b.^(j-1).*t.^(a*j-a).*sin(j.*pi.*rho);
end;
end

b2. avg_f1.txt:

function y=avg_f1(t,m,rho,r,a,b)

y=0;
x=r*sin(pi*rho);

b3. f2.txt:

function y=f2(t,m,rho,r,a,b,tr,nrho)
if nargin==7
nrho=size(tr,1)+1;
end;
if(rho<=0||rho>=1||b<=0||r<=0||a<=0||a>=2||r>=tr(rho*nrho,m))
y=1e7;
else
y=-sum(log(avg_f1(t,m,rho,r,a,b)));
end;
end

b4. main7.txt:

clear all;clc;
finp1='tr.xls';   finp2='AMEX.xls';
tr=xlsread(finp1,1,'B2:K10');
tr(tr>=1000)=10;
dat= xlsread(finp2,1,'D1:D1200');
[nr_tr nc_tr]=size(tr);
m=2;
tic;
for i=1:nr_tr
rho=i/(nr_tr+1);
lf=@(th) f2(dat,m,rho,th(1),th(2),th(3),tr);
[x fx]= patternsearch(lf,[min(tr(i,m)/2,.6),.4,.4],[],[],[],[],[0,0,0],[tr
(rho*10,m),2,10000000000])
    A(i,:)=[x,rho,fx];
fprintf('i = %3d\t,\ttime : %4d seconds\n',i,round(toc));
end;
[lmin,imin]=min(A(:,5));
th_hat=A(imin,:);
rh=th_hat(1);
ah=th_hat(2);
bh=th_hat(3);
rhoh=th_hat(4);
th_hat