Am I doing anything wrong when trying to calculate the expectation value in Qiskit on the real hardware?

I am trying to use the method of calculating the expectation values of Pauli operators in Qiskit which I found here.

However, the results obtained via IBMQ sampling differ significantly from the exact value, even though I perform measurement error mitigation.

The state, in which I am trying to calculate the expectation values, is prepared using a simple circuit

          ┌─────────────┐                   ┌───┐                     
q_0: ─────┤ RY(-1.8018) ├──■────────■───────┤ X ├────────────────────
     ┌───┐└──────┬──────┘┌─┴─┐      │       └─┬─┘                    
q_1: ┤ X ├───────■───────┤ X ├──────┼─────────┼──────────────────────
     └───┘               └───┘┌─────┴──────┐  │                 ┌───┐
q_2: ─────────────────────────┤ RY(2.2489) ├──■─────────■───────┤ X ├
                              └────────────┘     ┌──────┴──────┐└─┬─┘
q_3: ────────────────────────────────────────────┤ RY(0.99778) ├──■──
                                                 └─────────────┘     

which after transpiling looks as follows:

For reference, for each Pauli operator I print the exact expectation value and also calculate it a few times using the QASM simulator.

Here’s the code:

circuit = QuantumCircuit(4)
circuit.x(1)
circuit.cry(-1.80184863, 1, 0)
circuit.cx(0,1)
circuit.cry(2.24892942,0,2)
circuit.cx(2,0)
circuit.cry(0.9977846,2,3)
circuit.cx(3,2)
psi = CircuitStateFn( circuit )
    
paulis = [ Pauli([1,1,0,0],[1,1,0,0]),  Pauli([1,1,1,1],[1,0,0,1]) ]
shots = 8000
reps = 3   

backend_qasm = qiskit.Aer.get_backend( 'qasm_simulator' )
q_instance_qasm = QuantumInstance( backend_qasm, shots = shots )
    
load_account()
provider = get_provider( hub='ibm-q' )
backend_ibmq = least_busy( provider.backends(filters=lambda x: x.configuration().n_qubits >= 4 and not x.configuration().simulator) )
q_instance_ibmq = QuantumInstance( backend = backend_ibmq,
                                   shots = shots,
                                   measurement_error_mitigation_cls = CompleteMeasFitter,
                                   measurement_error_mitigation_shots = shots )
print(f'IBMQ backend: {backend_ibmq}.n')     

for pauli in paulis:
    print(f'Pauli operator: {pauli}.')
    pauli = WeightedPauliOperator([[1., pauli]]).to_opflow()

    measurable_expression = StateFn( pauli, is_measurement = True ).compose( psi )
    expectation = PauliExpectation().convert( measurable_expression )

    expect_exact = psi.adjoint().compose( pauli ).compose( psi ).eval().real

    print( f'Exact expectation value: {expect_exact}.' )
    for r in range(reps):
        sampler_qasm = CircuitSampler( q_instance_qasm ).convert( expectation )
        expect_sampling_qasm = sampler_qasm.eval().real
        print( f'Exact expectation, QASM sampling: {expect_sampling_qasm}.' )
    for r in range( reps ):
        sampler_ibmq = CircuitSampler( q_instance_ibmq ).convert( expectation )
        expect_sampling_ibmq = sampler_ibmq.eval().real
        print( f'Exact expectation, IBMQ sampling: {expect_sampling_ibmq}.' )

    print()

And here’s the output:

IBMQ backend: ibmq_ourense.

Pauli operator: IIYY.
WARNING - The skip Qobj validation does not work for IBMQ provider. Disable it.
Exact expectation value: -0.4201884924852.
Exact expectation, QASM sampling: -0.42275.
Exact expectation, QASM sampling: -0.4175.
Exact expectation, QASM sampling: -0.4165.
Exact expectation, IBMQ sampling: -0.19053720838213.
Exact expectation, IBMQ sampling: -0.33771371840093.
Exact expectation, IBMQ sampling: 0.14870401826006.

Pauli operator: YZZY.
Exact expectation value: 0.22895884311365.
Exact expectation, QASM sampling: 0.237.
Exact expectation, QASM sampling: 0.2385.
Exact expectation, QASM sampling: 0.2345.
Exact expectation, IBMQ sampling: 0.06862734682344.
Exact expectation, IBMQ sampling: 0.10246703115813.
Exact expectation, IBMQ sampling: 0.13078427261863.

Am I doing something conceptually wrong?

Is there an obvious way to improve the results? (except for doing more shots) Or is it what I should expect to get, given the device’s gate fidelities?

Any thoughts/suggestions/corrections greatly appreciated.