Digging deeper into the FitDist function I am using the script allfitdist to find which distribution best fits my data according to the tests included (BIC, AIK etc...) - these all stem from the script using the fitdist function in the statistics toolbox.
I manage to do this and get a nice historgram with the various distributions superimposed on my data. This question poses something similar to what I am struggling with, but I want to go a step further and know if it is possible to get the further parameters (for example the skew and kurtosis) of the output distributions.
I have literally no idea how I can do this. The output for a norm dist only contains mean and st.dev, but it seems from the figure that there is skewness and kurtosis - I wouldn't believe them to both be exactly zero either, given my data is market return of a stock (unless the fitdist function used has that as a assumption??)
If it would be helpful, I could post a picture of the output.
 A: Based on the documentation of fitdist the distribution argument used determines the class type of the returned probability distribution object. In the case of a normal distribution $N(\mu,\sigma^2)$, that means returning "only two" parameters, $\mu$ and $\sigma$, automatically (if you were using a $t$-distribution for example you expect getting the degrees of freedom $\nu$, etc...). So :
rng(123,'twister'); %Fix the seed
sample = randn(10^5,1)*2;
pdobj = fitdist(sample,'Normal')

pdobj = 

  NormalDistribution

  Normal distribution
       mu = 9.0349e-05   [-0.006102, 0.00628269]
    sigma =   0.999144   [0.994784, 1.00354]

As you can see yourself (and read in the documentation) $\mu$ and $\sigma$ are given as ML estimates so you could get them directly using mean and std. In the same manner, you can go ahead and use skewness and kurtosis directly to get the skewness and kurtosis of your data respectively. As such and going back to your original question: Yeap, it is possible to get further paramters. 
So for example assuming you are using the previous sample still, calculate the structure produced by allfitdist:
afobj = allfitdist(sample,'PDF');

then go ahead and find the index of afobj where the distribution name 'normal' is stored:
Indx = find(strcmp({afobj.DistName},'normal'));

add in that distribution parameter names two fields for skewness and kurtosis respectively:
afobj(Indx).ParamNames{3} = 'skewness';
afobj(Indx).ParamNames{4} = 'kurtosis';

and then you finally append the values of the skewness and the kurtosis in the field holding the parameter values:
afobj(Indx).Params(3:4) = [ skewness( sample), kurtosis(sample)];

Clearly you can do this hocus-pocus within allfitdist and not have to do that more that once ever again. So simply add:
Indx = find(strcmp({D.DistName},'normal'));
D(Indx).ParamNames{3} = 'skewness';
D(Indx).ParamNames{4} = 'kurtosis';
D(Indx).Params(3:4) = [ skewness( data), kurtosis(data)];

somewhere towards the end of the main function in the function file fidist.m (D is the structure's name within the function allfitdist). Clearly you can also add descriptions like:
afobj(Indx).ParamDescription{3} = '3rd moment'; 
afobj(Indx).ParamDescription{4} = '4th moment';

etc. etc. but that's a issue of personal tastes. A final thing that I almost forgot: with fitdist you can do the same qualitatives changes but it will be more cumbersome to actually implement them as one will need to go ahead and edit MATLAB class objects directly as the dot notation used in the previous example is not applicable there (strictly speaking, ParameterNames in a fitdist object are read-only properties and one cannot set them directly; thankfully the allfitdist result structure does not have these characteristics). 
