#include <iostream>
#include <bits/stdc++.h>
#define mp make_pair
#define pb push_back
using namespace std;
using ll = long long int;
const int p = 1e9 + 7;
ll modpower(ll base, ll power, ll mod=p){
    ll ans =1;
    base%=p;
    while(power){
        if(power&1){
            ans*=base;
            ans%=mod;
        }
        base*=base;
        base%=mod;
        power>>=1;
    }
    return ans;
}
template<typename T>
ostream& operator+(ostream& out, const vector<T> &vec){
    for(const auto &x : vec){
        out<<x<<" ";
    }
    out<<"\n";
    return out;
}
template<typename T>
ostream& operator*(ostream& out, const vector<T> &vec){
    for(const auto &x : vec){
        out+x;
    }
    return out;
}
template<typename T>
istream& operator>>(istream& in, vector<T> &vec){
    for(auto &x : vec){
        in>>x;
    }
    return in;
}
const int lgn = 20;
void solve(){
    int n;
    cin>>n;
    ++n;
    vector<vector<int>> adj(n);
    for(int i=1;i<n;i++){
        int p;
        cin>>p;
        adj[i].pb(p);
        adj[p].pb(i);
    }
    vector<vector<int>> queries(n);
    vector<array<int, lgn>> lift(n);
    vector<int> intime(n);
    vector<int> outtime(n);
    function<void(int,int,int, int&)> dfs = [&](int u,int par, int dep, int &t){
        lift[u][0] = par;
        intime[u] = t++;
        for(int i=1;i<lgn;i++){
            lift[u][i] = lift[lift[u][i-1]][i-1];
        }
        queries[dep].pb(u);
        for(const auto &v : adj[u]){
            if(v==par) continue;
            dfs(v, u, dep + 1, t);
        }
        outtime[u] = t++;
    };
    int time=0;
    dfs(0,0,0,time);
    function<bool(int,int)> ispar = [&](int p, int u){
        return intime[p]<=intime[u] && outtime[p]>=outtime[u];
    };
    function<int(int,int)> lca = [&](int u,int v){
        if(ispar(u,v)) return u;
        if(ispar(v,u)) return v;
        for(int i=lgn-1;i>=0;--i){
            if(!ispar(lift[u][i], v)){
                u = lift[u][i];
            }
        }
        return lift[u][0];
    };
    ll ans = 0;
    for(auto &q : queries){
        if((int)q.size()==0) break;
        if(q[0]==0){
            ans+=modpower(2, n-1);
            continue;
        }
        int y = q.size();
        for(int i=1;i<y;i++){
            q.pb(lca(q[i], q[i-1]));
        }
        q.pb(0);
        sort(q.begin(), q.end(), [&](const int &x,const int &y){return intime[x] < intime[y];});
        q.resize(unique(q.begin(), q.end()) - q.begin());
        int x = q.size();
        vector<vector<int>> vadj(x);
        vector<pair<int,int>> currtree = {{0,0}};
        for(int i=1;i<x;i++){
            while(!ispar(currtree.back().first, q[i])){
                currtree.pop_back();
            }
            vadj[currtree.back().second].pb(i);
            vadj[i].pb(currtree.back().second);
            currtree.pb({q[i], i});
        }
        vector<ll> dp0(x, 1);
        vector<ll> dp1(x, 1);
        function<void(int,int)> vdfs = [&](int u,int par){
            ll total = 1;
            ll none = 1;
            ll ratiosum = 0;
            if(vadj[u].size()==1 && vadj[u][0]==par) return;
            for(const auto &v : vadj[u]){
                if(v==par) continue;
                vdfs(v, u);
                total*=(dp0[v] + dp1[v]);
                total%=p;
                none*=dp0[v];
                none%=p;
                ratiosum+=dp1[v]*modpower(dp0[v], p-2);
                ratiosum%=p;
            }
            dp1[u] = none*ratiosum;
            dp1[u]%=p;
            dp0[u] = total - dp1[u];
            dp0[u]%=p;
            if(dp0[u] < 0) dp0[u]+=p;
        };
        vdfs(0, 0);
        ans+=dp1[0]*modpower(2, n - y);
        ans%=p;
    }
    cout<<ans<<"\n";
}
int main() {
    ios_base::sync_with_stdio(false);
    cin.tie(NULL);
    solve();
}

Measured memory_bytes: 60936