function [moments,pricestat,pctls,freqs]=v25_calcstatistics(STAT,XSIM,pctls_values,freqs_range)
%Function provides the following outpus:

%moments: MxNxOx8 array with M=Number of scenarios; N=Number of y-states (y-variation); O=Number of variables in model; 8=Number of moments calculated
%The moments are the following:Mean  Std. Dev. Skewness  Kurtosis  Min
%Max       %LB       %UB

%pricestat: MxNx5 array with M and N as before; 5=moments (min, mean-sd,mean, mean+sd, max)

%freqs: MxNxQxS with M and N as before; Q=Number of response variables;
%S=Number of frequency ranges to calculate; The freqs_range needs to be a
%nx2 matrix with n as number of frequency ranges and 2 numbers per row
%indicate the range (lowest and highest value).

%pctls: MxNxQxP array with M, N, and Q as before; P=Number of percentiles calculated;
%Inputs:
%XSIM(M).x(sim1,Q,sim2,N): M and N as before, sim1 and sim2 are simulation variables (dimensions)
moments=NaN(size(STAT,1),size(STAT,2),size(STAT(1).moments(:,1),1),8);
pricestat=NaN(size(STAT,1),size(STAT,2),5);
freqs=NaN(size(STAT,1),size(STAT,2),size(STAT(1).moments(:,1),1)-2,size(freqs_range,1));
pctls=NaN(size(STAT,1),size(STAT,2),size(STAT(1).moments(:,1),1)-2,size(pctls_values,2));
for i=1:size(STAT,1)
    tempdata=permute(reshape(permute(XSIM(i).x,[4 2 1 3]),size(XSIM(i).x,4),size(XSIM(i).x,2),[]),[3 2 1]);
    pctls(i,:,:,:)=permute(prctile(tempdata,pctls_values),[3 2 1]);
%       t1=permute(prctile(tempdata(:,:,:),pctls_values),[3 1 2]);size(t1)  

    for i2=1:size(freqs_range,1)
        freqs(i,:,:,i2)=permute(sum(tempdata >= freqs_range(i2,1) & tempdata <= freqs_range(i2,2) ),[3 2 1])/(size(XSIM(i).x,1)*size(XSIM(i).x,3));
    end;
%     tic;%Alternative way to compute this:
%     t1=arrayfun( @(x)sum(tempdata>=x-tol & tempdata<=x+tol)/size(XSIM(i).x,1)*size(XSIM(i).x,3), freqs_values,'UniformOutput', false);
%     toc;
    clear tempdata;
    
    for i2=1:size(STAT,2)
        moments(i,i2,:,:)=STAT(i,i2).moments(:,:);

        pricestat(i,i2,1)=STAT(i,i2).moments(9,5);
        pricestat(i,i2,2)=STAT(i,i2).moments(9,1)-STAT(i,i2).moments(9,2);
        pricestat(i,i2,3)=STAT(i,i2).moments(9,1);
        pricestat(i,i2,4)=STAT(i,i2).moments(9,1)+STAT(i,i2).moments(9,2);
        pricestat(i,i2,5)=STAT(i,i2).moments(9,6);
        
%         freqs(i,i2,:,:)=sum(XSIM(i).x(:,:,:,i2) == freqs_values)/numel(XSIM(i).x(:,:,:,i2));  
    end;
end;