use std::io::*;
use std::str::FromStr;
use std::collections::BinaryHeap;
use std::collections::BTreeMap;

#[allow(dead_code)]
fn read<T: FromStr>() -> T {
    let stdin = stdin();
    let stdin = stdin.lock();
    let token: String = stdin
        .bytes()
        .map(|c| c.expect("failed to read char") as char)
        .skip_while(|c| c.is_whitespace())
        .take_while(|c| !c.is_whitespace())
        .collect();
    token.parse().ok().expect("failed to parse token")
}

#[allow(dead_code)]
fn pow_mod(base: i64, power: i64, modulo: i64) -> i64 {
    let mut base = base;
    let mut power = power;
    let mut ans = 1;
    while power > 0 {
        if power & 1 == 1 { ans = ans * base % modulo; }
        power >>= 1;
        base = base * base % modulo;
    }
    ans
}

#[allow(dead_code)]
fn inv_mod(a: i64, modulo: i64) -> i64 {
    pow_mod(a, modulo - 2, modulo)
}

type Weight = i32;
#[allow(dead_code)]
#[derive(Clone, Debug, Eq, PartialEq)]
struct Edge {
    index: usize,
    src: usize,
    dest: usize,
    weight: Weight,
}
#[allow(dead_code)]
impl Edge {
    fn new(index: usize, src: usize, dest: usize, weight: Weight) -> Edge {
        Edge {
            index: index,
            src: src,
            dest: dest,
            weight: weight,
        }
    }
}
impl std::cmp::Ord for Edge {
    fn cmp(&self, other: &Self) -> std::cmp::Ordering {
        (self.weight, self.src, self.dest).cmp(&(other.weight, other.src, other.dest))
    }
}
impl std::cmp::PartialOrd for Edge {
    fn partial_cmp(&self, other: &Self) -> Option<std::cmp::Ordering> {
        Some(self.cmp(other))
    }
}
type Edges = Vec<Edge>;
type Graph = Vec<Edges>;
#[allow(dead_code)]
fn graph_new(n: usize) -> Graph {
    vec![vec![]; n]
}
#[allow(dead_code)]
fn add_bi_edge(graph: &mut Graph, index: &mut usize, src: usize, dest: usize, weight: Weight) {
    graph[src].push(Edge::new(*index, src, dest, weight));
    graph[dest].push(Edge::new(*index, dest, src, weight));
    *index += 1;
}
#[allow(dead_code)]
fn add_uni_edge(graph: &mut Graph, index: &mut usize, src: usize, dest: usize, weight: Weight) {
    graph[src].push(Edge::new(*index, src, dest, weight));
    *index += 1;
}


#[allow(dead_code)]
struct LCA {
    log_size: usize,
    parent: Vec<Vec<usize>>,
    depth: Vec<usize>,
}
#[allow(dead_code)]
impl LCA {
    fn new(g: &Graph, root: usize) -> LCA {
        let n = g.len();
        let mut log_size = 0;
        while (1 << log_size) < n { log_size += 1; }
        let mut parent = vec![vec![root; log_size]; n];
        let mut depth = vec![0; n];

        let mut que = vec![];
        que.push(root);
        while let Some(node) = que.pop() {
            for edge in g[node].iter() {
                if edge.dest == parent[node][0] { continue; }
                parent[edge.dest][0] = node;
                depth[edge.dest] = depth[node] + 1;
                que.push(edge.dest);
            }
        }

        for iter in 0..log_size - 1 {
            for i in 0..n {
                parent[i][iter + 1] = parent[parent[i][iter]][iter];
            }
        }
        LCA {
            log_size: log_size,
            parent: parent,
            depth: depth,
        }
    }

    fn lca(&self, mut u: usize, mut v: usize) -> usize {
        if self.depth[u] > self.depth[v] { std::mem::swap(&mut u, &mut v); }
        for i in 0..self.log_size {
            if (((self.depth[u] ^ self.depth[v]) >> i) & 1) == 1 { v = self.parent[v][i]; }
        }
        if u == v { return u; }
        for i in (0..self.log_size).rev() {
            if self.parent[u][i] != self.parent[v][i] {
                u = self.parent[u][i];
                v = self.parent[v][i];
            }
        }
        return self.parent[u][0];
    }
}

const MOD: i64 = 1e+9 as i64 + 7;
struct Solver {
    n: usize,
    g: Graph,
    lca: LCA,
}

impl Solver {
    fn new() -> Solver {
        let n: usize = read::<usize>() + 1;
        let mut g = graph_new(n);
        let mut e = 0;
        for i in 0..n-1 {
            let p: usize = read();
            add_bi_edge(&mut g, &mut e, p, i + 1, 1);
        }
        let lca = LCA::new(&g, 0);
        Solver {
            n: n,
            g: g,
            lca: lca,
        }
    }
    fn dfs(&self, from: usize, parent: usize, d: usize, depths: &mut Vec<Vec<usize>>) {
        if depths.len() == d { depths.push(vec![]); }
        depths[d].push(from);
        for edge in self.g[from].iter() {
            if edge.dest == parent { continue; }
            self.dfs(edge.dest, from, d + 1, depths);
        }
    }
    fn get_parent(&self, node: usize, parents: &mut BTreeMap<usize, usize>) -> usize {
        let mut ret = node;
        if parents.contains_key(&node) {
            let pnode = *parents.get(&node).unwrap();
            ret = self.get_parent(pnode, parents);
            parents.insert(node, ret);
        }
        ret
    }
    fn merge(&self, childs: &Vec<usize>, cnts: &BTreeMap<usize, (i64, i64)>, parent: usize, parents: &mut BTreeMap<usize, usize>) -> (i64, i64) {
        let mut child_cnts = vec![];
        let mut all = 1;
        let mut zeros = 1;
        for child in childs {
            let p = &self.get_parent(*child, parents);
            if *p == parent { continue; }
            let cnt = cnts.get(p).unwrap();
            child_cnts.push(cnt);
            all = (all * (cnt.0 + cnt.1)) % MOD;
            zeros = (zeros * cnt.0) % MOD;
            parents.insert(*p, parent);
            // println!("{} {:?}", p, cnt);
        }
        let mut ones = 0;
        for cnt in child_cnts {
            ones = (ones + zeros * inv_mod(cnt.0, MOD) % MOD * cnt.1 % MOD) % MOD;
        }
        ((all - ones + MOD) % MOD, ones)
    }
    fn count_cnt(&self, nodes: &Vec<usize>) -> i64 {
        // println!("depth {:?}", nodes);
        let mut cnts = BTreeMap::<usize, (i64, i64)>::new();
        let mut childs = BTreeMap::<usize, Vec<usize>>::new();
        let mut que = BinaryHeap::<(usize, usize)>::new();
        for node in nodes { cnts.insert(*node, (1, 1)); }
        for pair in nodes.windows(2) {
            let pnode = self.lca.lca(pair[0], pair[1]);
            let pnode_childs = childs.entry(pnode).or_insert(vec![]);
            if pnode_childs.last().is_none() || *pnode_childs.last().unwrap() != pair[0] {
                pnode_childs.push(pair[0]);
            }
            pnode_childs.push(pair[1]);
            que.push((self.lca.depth[pnode], pnode));
        }
        let mut last_cnt = 1;
        let mut parents = BTreeMap::<usize, usize>::new();
        while let Some((_depth, node)) = que.pop() {
            if cnts.contains_key(&node) { continue; }
            let child = childs.get(&node).unwrap();
            // println!("cnts:   {:?}", cnts);
            // println!("parents {:?}", parents);
            // println!("childs  {:?}", child);
            let cnt = self.merge(child, &cnts, node, &mut parents);
            cnts.insert(node, cnt);
            // println!("{} {:?}", node, cnt);
            last_cnt = cnt.1
        }
        // println!("");
        last_cnt
    }
    fn solve(&mut self) {
        let mut pow2 = vec![1; self.n + 1];
        for i in 0..self.n {
            pow2[i + 1] = pow2[i] * 2 % MOD;
        }
        let mut depths = vec![];
        self.dfs(0, 0, 0, &mut depths);
        // println!("{:?}", depths);
        // println!("{:?}", self.lca.parent);
        // println!("{:?}", self.lca.depth);
        let mut ans = 0;
        for nodes in depths {
            let cnt = self.count_cnt(&nodes);
            let amp = pow2[self.n - nodes.len()];
            // println!("{} {}", cnt, amp);
            ans = (ans + cnt * amp) % MOD;
        }
        println!("{}", ans);
    }
}

fn main() {
    let mut solver = Solver::new();
    solver.solve();
}

Measured memory_bytes: 112892