function[nmb,nme,mb,me,rmse,corr] = calcStats(observed,predicted)

    nmb = [];
    nme = [];
    mb = [];
    me = [];
    rmse = [];
    corr = [];
    for vector = 1:length(observed)
        obs = observed{vector};
        predics = predicted{vector};

        % Removing NaN values
        obs_n = obs(~isnan(obs));
        predic_n = predics(~isnan(obs));
        %if including depolarization
        nans_left = isnan(predic_n);
        obs_n(nans_left) = [];
        predic_n(nans_left) = [];
        
        % Calculating NMB and NME
        nmb{vector} = (mean(predic_n) - mean(obs_n)) / mean(obs_n);
        nme{vector} = mean(abs(predic_n - obs_n)) / mean(obs_n);
        
        % Calculating MB and ME
        mb{vector} = mean(predic_n - obs_n);
        me{vector} = mean(abs(predic_n - obs_n));
        
        % Calculating RMSE
        rmse{vector} = sqrt(mean((predic_n - obs_n).^2));
        
        % Calculating CORR
        corr{vector} = corrcoef(predic_n, obs_n);
        corr{vector} = corr{vector}(1, 2);

    end

