diff --git a/referenceqvm/qvm_unitary.py b/referenceqvm/qvm_unitary.py
index 1114c47..fcd526e 100644
--- a/referenceqvm/qvm_unitary.py
+++ b/referenceqvm/qvm_unitary.py
@@ -21,7 +21,10 @@
 from pyquil.quil import Program
 from pyquil.quilbase import *
 
-from referenceqvm.unitary_generator import tensor_gates
+import numpy as np
+import scipy.sparse as sps
+
+from referenceqvm.unitary_generator import tensor_gates, lifted_gate
 from referenceqvm.qam import QAM
 
 
@@ -108,5 +111,33 @@ def expectation(self, pyquil_program, operator_programs=[Program()]):
         :return: expectation value of the operators.
         :rtype: float
         """
-        # TODO
-        raise NotImplementedError()
+        # get the largest amount of qubits manipulated in a program, so that, later,
+        # all the matrices can be lifted to process the full qubit system
+        all_programs = [pyquil_program] + operator_programs
+        highest_qubit_index = max({max(program.get_qubits() or {0}) for program in all_programs})
+        num_qubits = highest_qubit_index + 1
+        
+        # get the evolved state of the preparatory program
+        prep_umat = self.unitary(pyquil_program)
+        prep_umat = lifted_gate(0, prep_umat, num_qubits)
+
+        # the ground state vector
+        ground_state = sps.csc_matrix([1.0], shape=(2**num_qubits,1))
+
+
+        # get the preperatory state by multiplying the preperatory program's
+        # matrix against the ground state assumed by quil programs
+        prep_state = prep_umat.dot(ground_state)
+        conj_prep_state = prep_state.getH()
+
+        # obtain unitary form of operators
+        observables = [self.unitary(program) for program in operator_programs]
+        observables = [lifted_gate(0, umat, num_qubits) for umat in observables]
+
+        # calculate the expectation value of the operators given the prep state
+        results = [conj_prep_state.dot(observable.dot(prep_state)) for observable in observables]
+
+        # getting the float value of the 1x1 sparse matrices
+        results = [result.sum() for result in results]
+
+        return results
diff --git a/referenceqvm/tests/test_unitary.py b/referenceqvm/tests/test_unitary.py
index 5177bfe..8b2f40d 100644
--- a/referenceqvm/tests/test_unitary.py
+++ b/referenceqvm/tests/test_unitary.py
@@ -1,5 +1,6 @@
 from pyquil.quil import Program
 from pyquil.gates import *
+from pyquil.paulis import sZ, sI, sX
 import numpy as np
 import pytest
 from referenceqvm.gates import gate_matrix
@@ -68,3 +69,41 @@ def test_unitary_errors(qvm_unitary):
     prog.inst(("hello2", 0))
     with pytest.raises(TypeError):
         qvm_unitary.unitary(prog)
+
+def test_unitary_expectation_identity(qvm_unitary):
+    """
+    Tests that the unitary QVM calculates the expectation value of the
+    identity observable on a ground state correctly.
+    """
+    identity = Program()
+    pauli_operators = [sI(0)]
+
+    expects = qvm_unitary.expectation(identity, [observable.program for observable in pauli_operators])
+    
+    assert len(expects) == 1
+    assert np.allclose(expects, [1])
+
+
+def test_unitary_expectation_bell_state(qvm_unitary):
+    """
+    Tests that the unitary QVM calculates the expectation value of various pauli operators
+    on the bell state correctly
+    """
+    bell_state = Program(
+        H(0),
+        CNOT(0, 1),
+    )
+    pauli_operators = [
+        sZ(0) * sZ(1),
+        sZ(0),
+        sZ(1),
+        sX(0) * sX(1),
+    ]
+    expects = qvm_unitary.expectation(bell_state, [observable.program for observable in pauli_operators])
+
+    assert len(expects) == 4
+    assert np.allclose(expects, [1, 0, 0, 1])
+
+
+
+
diff --git a/requirements.txt b/requirements.txt
index 07be7f6..8f60011 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,6 +1,6 @@
 numpy>=1.11.1
 scipy>=0.18.1
-pyquil>=1.4.2
+pyquil==1.9.0
 
 # For testing
 pytest>= 3.0.0