Skewness of fitted mixture not correct? I fitted a gaussian mixture to my financial data. 
The values are:
$\pi= 0.3$
$\mu_1= -0.01$
$\mu_2= 0.01$
$\sigma_1=0.01$
$\sigma_2=0.03$
One can see, that both single distributions have a mean of almost zero, wherease one has a high volatility and the other a low volatility. The normal distribution 1, the green one with the high peak has the parameters $\mu_1$ and $\sigma_1$ and occurs (this is pi from output of normalmixEM) with a probability of 0.39. The normal distribution 2 with the smaller peak and the higher volatility has the parameters $\mu_2$ and $\sigma_2$ and a probability of $1-0.39$.
I imagine the generating of the mixture density as follows:
We have a distribution which is quite probable ($\pi=1-0.39)$ and has $\mu_2$., if the mixture density is done, we "add" a second distribution which is a bit shifted to the left (this one occurs with a probability of 0.39 and has a negative mean). Since the distribution we add lies a bit more to the left I would expect, that the mixture density has a negative skew, since the left tail of the resulting mixture will be heavier?
I control this, which gives a positive skew of 0.7. Now my question is: Why? I would expect a negative skew, since I thought the mixture density will have a a fatter left tail, since we add to the probable distribution with positive mean a second distribution which is a bit shifted to the left?
 A: I think it is time to start paying attention to the answers you receive:  https://stats.stackexchange.com/a/54734

It has been largely discussed the validity of using a measure of skewness in multimodal distributions, since its interpretation becomes unclear. This is the case of finite mixtures. If your mixture looks (or is) unimodal, then you can use this value to understand a bit how asymmetric it is.

In addition, note that MLE is based on maximising the likelihood, not on matching moments.
Check the following example using a simulated sample from a skew-normal (the sample skewness and the skewness of the mixture are different). 
library(sn)
library(moments)
library(mixtools)

set.seed(4321)    
samp <- rsn(1000,0,1,10)

skewness(samp)

mix<-normalmixEM(samp,k=2,fast=TRUE)

# Sampling from a 2-gaussian mixture
gaussmix <- function(n,m1,m2,s1,s2,alpha) {
    I <- runif(n)<alpha
    rnorm(n,mean=ifelse(I,m1,m2),sd=ifelse(I,s1,s2))
}

# A simulated sample
samp2 <- gaussmix(100000,mix$mu[1],mix$mu[2],mix$sigma[1],mix$sigma[2],mix$lambda[1])

skewness(samp2)

Update
The measure of skewness based on the third moment is driven by the tails, rather than the level of asymmetry of the distribution. For example, a Student-$t$ distribution with $\nu<3$ degrees of freedom has undefined skewness in spite of being clearly symmetric. For this reason, it is preferred to use a quantile-based measure of skewness. An example of this of this is the AG measure of skewness which can be defined for any unimodal distribution $F$ as
$$AG=1-2F(\mbox{mode}).$$
This can be estimated by using a nonparametric estimator of $F$, such as a kernel estimator as follows
rm(list=ls())
library(sn)
library(moments)
library(mixtools)

# AG measure of skewness: insert the data and the interval c(minim,maxim) where the mode is located

AG <- function(data,minim,maxim){
n = length(data)
hb = (4*sqrt(var(data))^5/(3*n))^(1/5)

     kn = function(x){
     k = r = length(x)
     for(i in 1:k) r[i] = mean(dnorm((x[i]-data)/hb))/hb
     return(r)
      } 

    mode = optimise(f = kn, interval = c(minim,maxim),maximum = TRUE)$maximum

     KN = function(x){
     k = r = length(x)
    for(i in 1:k) r[i] = mean(pnorm((x[i]-data)/hb))
    return(r)
    } 
    return(1-2*KN(mode))
}


# Simulated sample
set.seed(707)
samp <- rsn(1000,0,1,10)

# AG and moment-based skewness measure of the sample
AG(samp,0,2)
skewness(samp)

mix<-normalmixEM(samp,k=2,fast=TRUE)

# Sampling from a 2-gaussian mixture
gaussmix <- function(n,m1,m2,s1,s2,alpha) {
    I <- runif(n)<alpha
    rnorm(n,mean=ifelse(I,m1,m2),sd=ifelse(I,s1,s2))
}

# A simulated sample
mix<-normalmixEM(samp,k=2,fast=TRUE)
samp2 <- gaussmix(100000,mix$mu[1],mix$mu[2],mix$sigma[1],mix$sigma[2],mix$lambda[1])

# AG and moment-based skewness measure of the fitted mixture
AG(samp2,0,2)
skewness(samp2)

As you can see, the AG measure of skewness is similar for both distributions, as expected.
The moral of the story is: The moment-based measure of skewness is NOT the best choice.
The same happens if this measure is applied to your data:
library(sn)
library(moments)
library(mixtools)

# AG measure of skewness: insert the data and the interval c(minim,maxim) where the mode is located

AG <- function(data,minim,maxim){
n = length(data)
hb = (4*sqrt(var(data))^5/(3*n))^(1/5)

     kn = function(x){
     k = r = length(x)
     for(i in 1:k) r[i] = mean(dnorm((x[i]-data)/hb))/hb
     return(r)
      } 

    mode = optimise(f = kn, interval = c(minim,maxim),maximum = TRUE)$maximum

     KN = function(x){
     k = r = length(x)
    for(i in 1:k) r[i] = mean(pnorm((x[i]-data)/hb))
    return(r)
    } 
    return(1-2*KN(mode))
}



# AG and moment-based skewness measure of the sample
AG(dat,-.1,.1)
skewness(dat)

# Sampling from a 2-gaussian mixture
gaussmix <- function(n,m1,m2,s1,s2,alpha) {
    I <- runif(n)<alpha
    rnorm(n,mean=ifelse(I,m1,m2),sd=ifelse(I,s1,s2))
}

# A simulated sample
mix<-normalmixEM(dat,k=2,fast=TRUE)
samp2 <- gaussmix(100000,mix$mu[1],mix$mu[2],mix$sigma[1],mix$sigma[2],mix$lambda[1])

# AG and moment-based skewness measure of the fitted mixture
AG(samp2,-.1,.1)
skewness(samp2)

